TcRay21qb64v2.txt

'TCRAY21 by Jark 2003
'Ported to QB64PE by Antoni Gual 2026




' Shared variables for SVGA management


Dim Shared xPlotMin%, xPlotMax%
Dim Shared yPlotMin%, yPlotMax%



Common Shared SizeScreen As Integer
Common Shared ScreenRes$, Name$
Common Shared ScrWidth As Long, ScrHeight As Long
Common Shared Red%, Green%, Blue%
Common Shared xScrCenter%, yScrCenter%

Const eps = 9.98012604599318D-322 ' See Antoni's version
'CONST eps = 1.401298E-45             ' Avoids logarithm overflow

' Constants for the TC-Lib programs
Const SaveDir$ = "" ' For files management
Const Pi = 3.14159265358979# ' For angles management
Const SQR12 = .7071068 'SQR(1/2)
Const SQR16 = .4082483 'SQR(1/6)
Const SQR23 = .8164966 'SQR(2/3)
Const SatCoeff = 84.51663 '100/ATN(SQR(6))
Const SatCoef = 1.183199E-02 '1/SatCoeff
Const PiOver180 = .0174533
Const PiOver2 = 1.5707963267949#

Const Gamma = 1.4, Aperture = 1



Dim Shared xCam(3), xLight(3), xCenter(3), UScreen(3), VScreen(3)

Type shapetyp
    Degree As Integer
    ShapeType As String * 2
    NumType As Integer
    HueRef As Single
    SatRef As Single
    x0 As Single
    y0 As Single
    z0 As Single
    dRef1 As Single
    dRef2 As Single
    Brilliance As Single
    Reflect As Single
    Texture As Integer
    Bump As Integer
End Type

Dim Shared NbObjects%
Dim Shared shape(64) As shapetyp, Object(64, 35), Matrix(64, 3, 3)
Dim Shared cReflect, cShadow, AntiAlias
Dim Shared t!


' Shared variables for the Perlin noises
Const octaves% = 7
Dim Shared Seed&, Noise%(129, 129), Persistence!, Span%, Amplitude!(16)
' Powers of two precalculation
Dim Shared Power2%(15)
For k% = 0 To 14: Power2%(k%) = 2 ^ k%: Next k%: Power2%(15) = &H8000

Clear
SizeScreen = 2

top:
NbObjects% = 0
Erase shape

ScrWidth = 1024
ScrHeight = 768
Screen _NewImage(ScrWidth, ScrHeight, 32)
_ScreenMove _Middle
_Title "Jark's TCRay Raytracer"
xPlotMin% = 0
xPlotMax% = ScrWidth - 1
yPlotMin% = 0
yPlotMax% = ScrHeight - 1
xScrCenter% = Int((xPlotMin% + xPlotMax%) / 2 + .5)
yScrCenter% = Int((yPlotMin% + yPlotMax%) / 2 + .5)
_Font _LoadFont("cour.ttf", 24)



Menu

t! = Timer
Draw.Space
t! = Timer - t!
Sleep
Do While Len(InKey$): Loop

GoTo top:

EndProg


Function Acos (x#)
    ' Returns the ArcCosinus of a real number

    If x# = 0 Then Acos = PiOver2: Exit Function

    If x# > 0 Then
        Acos = Atn(Sqr(1 - x# * x#) / x#)
    Else
        Acos = Pi + Atn(Sqr(1 - x# * x#) / x#)
    End If
End Function

Function Arg (xa, ya)
    ' Returns the argument (in radians) of a point in the (x,y) plan
    ' Ranges from -Pi to Pi

    If xa = 0 And ya = 0 Then Arg = 0: Exit Function
    If xa = 0 And ya >= 0 Then Arg = Pi / 2: Exit Function
    If ya = 0 And xa < 0 Then Arg = Pi: Exit Function
    If xa = 0 And ya < 0 Then Arg = -Pi / 2: Exit Function
    If xa > 0 Then Arg = Atn(ya / xa): Exit Function
    If xa < 0 And ya >= 0 Then Arg = Pi - Atn(-ya / xa): Exit Function
    If xa < 0 And ya < 0 Then Arg = -Pi + Atn(ya / xa): Exit Function

End Function


Function CosD (A) Static
    ' Corrects the wrong values provided by thea standard COS function
    ' a is in degrees

    A = A - Int(A / 360 + .5) * 360

    If A >= 0 And A <= 90 Then
        x# = A / 180 * Pi
        If Cos(x#) < .5 Then
            CosD = Sin(Pi / 2 - x#)
        Else
            CosD = Cos(x#)
        End If
        Exit Function
    End If

    If A < 0 And A >= -90 Then
        x# = -A / 180 * Pi
        If Cos(x#) < .5 Then
            CosD = Sin(Pi / 2 - x#)
        Else
            CosD = Cos(x#)
        End If
        Exit Function
    End If

    If A > 90 And A <= 180 Then
        x# = (180 - A) / 180 * Pi
        If Cos(x#) < .5 Then
            CosD = Sin(x# - Pi / 2)
        Else
            CosD = -Cos(x#)
        End If
        Exit Function
    End If

    If A < -90 And A >= -180 Then
        x# = (A - 180) / 180 * Pi
        If Cos(x#) < .5 Then
            CosD = -Sin(Pi / 2 - x#)
        Else
            CosD = -Cos(x#)
        End If
        Exit Function
    End If

End Function

Sub EndProg

    ' Quit the program
    Cls
    Print
    Print
    Print
    Print
    Print "               ÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ»"
    Print "               º         Thanks for using the TC-Lib            º"
    Print "               º  Please report any bug, remark, comment... to  º"
    Print "               º          mandelbrot.dazibao@free.fr            º"
    Print "               º    and feel free to leave your footprint on    º"
    Print "               º       http://mandelbrot.dazibao.free.fr        º"
    Print "               ÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ"
    _Display
    Sleep
    End
End Sub



Sub FullScreen
    ' The Pset24 routine can plot any pixel on the screen

    xPlotMin% = 0
    xPlotMax% = ScrWidth - 1
    yPlotMin% = 0
    yPlotMax% = ScrHeight - 1
    xScrCenter% = Int((xPlotMin% + xPlotMax%) / 2 + .5)
    yScrCenter% = Int((yPlotMin% + yPlotMax%) / 2 + .5)
End Sub

Function Gamma.Correction (Value) Static
    ' Corrects the monitor rendering for dark pixels
    ' The Gamma correction coefficient should be set once and for all
    ' for you own computer monitor
    ' Gamma = 1 > no correction
    Gamma.Correction = ((Value / 100) ^ (1 / Gamma)) * 100
End Function

Sub HSVtoRGB (Hue, Sat, Value) Static
    ' ======================================================================
    ' Optimised version (35% faster than the previous version of the TC-Lib)
    ' ======================================================================

    ' Converts a HSV colour definition into RGB values via exact
    ' trigonometry calculations. Refer to HSVtoRGB1 for a simplified,
    ' faster but poorer conversion algorithm.
    ' Hue ranges from 0 to 360 degrees. Red corresponds to Hue = 0
    ' Saturation ranges from 0 to 100
    ' Value ranges from 0 to 100
    ' Saturation = 0 provides grey tones, from black (Value = 0)
    ' to white (Value = 100).
    ' Red%, Green% and Blue% are shared variables
    ' The routine uses 4 constants:
    '    SQR(1/2) = .7071068
    '    SQR(1/6) = .4082483
    '    ATN(SQR(6)/100 = 1.183199E-02
    '    Pi/180 = .0174533

    Angle = (Hue - 150) * PiOver180
    Ur = Value * 2.55
    Radius = Ur * Tan(Sat * SatCoef)
    Vr = Radius * Cos(Angle) * SQR12
    Wr = Radius * Sin(Angle) * SQR16

    Red% = Ur - Vr - Wr
    Green% = Ur + Vr - Wr
    Blue% = Ur + Wr + Wr

    If Red% < 0 Then
        Rdim = Ur / (Vr + Wr)
        Red% = 0
        Green% = Ur + (Vr - Wr) * Rdim
        Blue% = Ur + 2 * Wr * Rdim
        GoTo Ctrl255
    End If

    If Green% < 0 Then
        Rdim = -Ur / (Vr - Wr)
        Red% = Ur - (Vr + Wr) * Rdim
        Green% = 0
        Blue% = Ur + 2 * Wr * Rdim
        GoTo Ctrl255
    End If

    If Blue% < 0 Then
        Rdim = -Ur / (Wr + Wr)
        Red% = Ur - (Vr + Wr) * Rdim
        Green% = Ur + (Vr - Wr) * Rdim
        Blue% = 0
        GoTo Ctrl255
    End If

    Ctrl255:
    If Red% > 255 Then
        Rdim = (Ur - 255) / (Vr + Wr)
        Red% = 255
        Green% = Ur + (Vr - Wr) * Rdim
        Blue% = Ur + 2 * Wr * Rdim
    End If

    If Green% > 255 Then
        Rdim = (255 - Ur) / (Vr - Wr)
        Red% = Ur - (Vr + Wr) * Rdim
        Green% = 255
        Blue% = Ur + 2 * Wr * Rdim
    End If

    If Blue% > 255 Then
        Rdim = (255 - Ur) / (Wr + Wr)
        Red% = Ur - (Vr + Wr) * Rdim
        Green% = Ur + (Vr - Wr) * Rdim
        Blue% = 255
    End If

End Sub

Function Module (x, y, z)
    ' Returns the module of a 3D vector
    Module = Sqr(x * x + y * y + z * z)
End Function

Sub RGBtoHSV (Hue, Saturation, Value) Static
    ' Returns the HSV codes of the current RGB combination
    Temp = (Red% + Green% + Blue%) / 3
    xa = (Green% - Red%) * SQR12
    ya = (Blue% + Blue% - Red% - Green%) * SQR16
    Hue = Arg(xa, ya) / PiOver180 + 150
    Saturation = Arg(Temp, Module(Red% - Temp, Green% - Temp, Blue% - Temp)) * SatCoeff
    Value = Temp / 2.55

    If Saturation = 0 Or Value = 0 Then Hue = 0
    If Hue < 0 Then Hue = Hue + 360
    If Hue >= 360 Then Hue = Hue - 360
End Sub

Function ScalProd (u(), v())
    ' Returns the scalar product of two vectors
    ScalProd = u(1) * v(1) + u(2) * v(2) + u(3) * v(3)
End Function


Function SinD (A)
    ' Correct Sine value for an angle in degrees
    SinD = CosD(90 - A)
End Function

Function Sqr3 (x#)
    ' Cubic root of a double length real number
    If x# = 0 Then Sqr3 = 0: Exit Function
    If x# > 0 Then
        Sqr3 = Exp(Log(x#) / 3)
    Else
        Sqr3 = -Exp(Log(-x#) / 3)
    End If
End Function

'




' ==============================================================

Sub Calc.Normal (Id%, xh!(), xn!(), Ray!()) Static
    ' Get the normal vector xn on object number Id%
    ' The hit point is xh is the object axis system
    ' The normal vector is returned in the main axis system for brilliance calculation
    ' xo is the same vector, but in the object axis system
    ' The ray is defined in the main axis system

    Dim xo!(3)

    '=======================================================
    ' Calculate the normal vector in object axis system

    Select Case shape(Id%).Degree
        Case 1
            xo!(1) = Object(Id%, 32)
            xo!(2) = Object(Id%, 33)
            xo!(3) = Object(Id%, 34)

        Case 2
            xo!(1) = (2 * Object(Id%, 26) * xh!(1) + Object(Id%, 29) * xh!(2) + Object(Id%, 31) * xh!(3) + Object(Id%, 32))
            xo!(2) = (2 * Object(Id%, 27) * xh!(2) + Object(Id%, 29) * xh!(1) + Object(Id%, 30) * xh!(3) + Object(Id%, 33))
            xo!(3) = (2 * Object(Id%, 28) * xh!(3) + Object(Id%, 30) * xh!(2) + Object(Id%, 31) * xh!(3) + Object(Id%, 34))

        Case 3
            xo!(1) = 3 * Object(Id%, 17) * xh!(1) ^ 2 + Object(Id%, 16) * xh!(2) * xh!(3) + 2 * Object(Id%, 20) * xh!(1) * xh!(2) + 2 * Object(Id%, 21) * xh!(1) * xh!(3) + Object(Id%, 23) * xh!(2) ^ 2 + Object(Id%, 24) * xh!(3) ^ 2
            xo!(1) = xo!(1) + 2 * Object(Id%, 26) * xh!(1) + Object(Id%, 29) * xh!(2) + Object(Id%, 31) * xh!(3) + Object(Id%, 32)
            xo!(2) = 3 * Object(Id%, 18) * xh!(2) ^ 2 + Object(Id%, 16) * xh!(3) * xh!(1) + 2 * Object(Id%, 22) * xh!(2) * xh!(3) + 2 * Object(Id%, 23) * xh!(2) * xh!(1) + Object(Id%, 25) * xh!(3) ^ 2 + Object(Id%, 20) * xh!(1) ^ 2
            xo!(2) = xo!(2) + 2 * Object(Id%, 27) * xh!(2) + Object(Id%, 30) * xh!(3) + KYX * xh!(1) + Object(Id%, 33)
            xo!(3) = 3 * Object(Id%, 19) * xh!(3) ^ 2 + Object(Id%, 16) * xh!(1) * xh!(2) + 2 * Object(Id%, 24) * xh!(3) * xh!(1) + 2 * Object(Id%, 25) * xh!(3) * xh!(2) + Object(Id%, 21) * xh!(1) ^ 2 + Object(Id%, 22) * xh!(2) ^ 2
            xo!(3) = xo!(3) + 2 * Object(Id%, 28) * xh!(3) + Object(Id%, 31) * xh!(1) + Object(Id%, 30) * xh!(2) + Object(Id%, 34)

        Case 4
            xo!(1) = 4 * Object(Id%, 1) * xh!(1) * xh!(1) * xh!(1)
            xo!(1) = xo!(1) + 2 * Object(Id%, 10) * xh!(1) * xh!(2) * xh!(2)
            xo!(1) = xo!(1) + 2 * Object(Id%, 11) * xh!(1) * xh!(3) * xh!(3)
            xo!(1) = xo!(1) + 3 * Object(Id%, 4) * xh!(1) * xh!(1) * xh!(2)
            xo!(1) = xo!(1) + 3 * Object(Id%, 5) * xh!(1) * xh!(1) * xh!(3)
            xo!(1) = xo!(1) + Object(Id%, 7) * xh!(2) * xh!(2) * xh!(2)
            xo!(1) = xo!(1) + Object(Id%, 8) * xh!(3) * xh!(3) * xh!(3)
            xo!(1) = xo!(1) + 3 * Object(Id%, 17) * xh!(1) * xh!(1)
            xo!(1) = xo!(1) + Object(Id%, 16) * xh!(2) * xh!(3)
            xo!(1) = xo!(1) + 2 * Object(Id%, 20) * xh!(1) * xh!(2)
            xo!(1) = xo!(1) + 2 * Object(Id%, 21) * xh!(1) * xh!(3)
            xo!(1) = xo!(1) + Object(Id%, 23) * xh!(2) * xh!(2)
            xo!(1) = xo!(1) + Object(Id%, 24) * xh!(3) * xh!(3)
            xo!(1) = xo!(1) + 2 * Object(Id%, 26) * xh!(1)
            xo!(1) = xo!(1) + Object(Id%, 29) * xh!(2)
            xo!(1) = xo!(1) + Object(Id%, 31) * xh!(3)
            xo!(1) = xo!(1) + Object(Id%, 32)

            xo!(2) = 4 * Object(Id%, 2) * xh!(2) * xh!(2) * xh!(2)
            xo!(2) = xo!(2) + 2 * Object(Id%, 12) * xh!(2) * xh!(3) * xh!(3)
            xo!(2) = xo!(2) + 2 * Object(Id%, 10) * xh!(2) * xh!(1) * xh!(1)
            xo!(2) = xo!(2) + 3 * Object(Id%, 6) * xh!(2) * xh!(2) * xh!(3)
            xo!(2) = xo!(2) + 3 * Object(Id%, 7) * xh!(2) * xh!(2) * xh!(1)
            xo!(2) = xo!(2) + Object(Id%, 9) * xh!(3) * xh!(3) * xh!(3)
            xo!(2) = xo!(2) + Object(Id%, 4) * xh!(1) * xh!(1) * xh!(1)
            xo!(2) = xo!(2) + 3 * Object(Id%, 18) * xh!(2) * xh!(2)
            xo!(2) = xo!(2) + Object(Id%, 16) * xh!(3) * xh!(1)
            xo!(2) = xo!(2) + 2 * Object(Id%, 22) * xh!(2) * xh!(3)
            xo!(2) = xo!(2) + 2 * Object(Id%, 23) * xh!(2) * xh!(1)
            xo!(2) = xo!(2) + Object(Id%, 25) * xh!(3) * xh!(3)
            xo!(2) = xo!(2) + Object(Id%, 20) * xh!(1) * xh!(1)
            xo!(2) = xo!(2) + 2 * Object(Id%, 27) * xh!(2)
            xo!(2) = xo!(2) + Object(Id%, 30) * xh!(3)
            xo!(2) = xo!(2) + KYX * xh!(1)
            xo!(2) = xo!(2) + Object(Id%, 33)

            xo!(3) = 4 * Object(Id%, 3) * xh!(3) * xh!(3) * xh!(3)
            xo!(3) = xo!(3) + 2 * Object(Id%, 11) * xh!(3) * xh!(1) * xh!(1)
            xo!(3) = xo!(3) + 2 * Object(Id%, 12) * xh!(3) * xh!(2) * xh!(2)
            xo!(3) = xo!(3) + 3 * Object(Id%, 8) * xh!(3) * xh!(3) * xh!(1)
            xo!(3) = xo!(3) + 3 * Object(Id%, 9) * xh!(3) * xh!(3) * xh!(2)
            xo!(3) = xo!(3) + Object(Id%, 5) * xh!(1) * xh!(1) * xh!(1)
            xo!(3) = xo!(3) + Object(Id%, 6) * xh!(2) * xh!(2) * xh!(2)
            xo!(3) = xo!(3) + 3 * Object(Id%, 19) * xh!(3) * xh!(3)
            xo!(3) = xo!(3) + Object(Id%, 16) * xh!(1) * xh!(2)
            xo!(3) = xo!(3) + 2 * Object(Id%, 24) * xh!(3) * xh!(1)
            xo!(3) = xo!(3) + 2 * Object(Id%, 25) * xh!(3) * xh!(2)
            xo!(3) = xo!(3) + Object(Id%, 21) * xh!(1) * xh!(1)
            xo!(3) = xo!(3) + Object(Id%, 22) * xh!(2) * xh!(2)
            xo!(3) = xo!(3) + 2 * Object(Id%, 28) * xh!(3)
            xo!(3) = xo!(3) + Object(Id%, 31) * xh!(1)
            xo!(3) = xo!(3) + Object(Id%, 30) * xh!(2)
            xo!(3) = xo!(3) + Object(Id%, 34)

        Case Else
    End Select

    ' Normalise the vector
    ModNorm = Module(xo!(1), xo!(2), xo!(3))
    If ModNorm = 0 Then Exit Sub
    For k% = 1 To 3: xo!(k%) = xo!(k%) / ModNorm: Next k%

    ' Convert to main axis system
    For k% = 1 To 3
        xn!(k%) = xo!(1) * Matrix(Id%, 1, k%) + xo!(2) * Matrix(Id%, 2, k%) + xo!(3) * Matrix(Id%, 3, k%)
    Next k%

    ' Reverse the vector if necessary
    Side = Ray!(4) * xn!(1) + Ray!(5) * xn!(2) + Ray!(6) * xn!(3)
    If Side < 0 Then
        For k% = 1 To 3: xn!(k%) = -xn!(k%): Next k%
    End If

End Sub

Sub Draw.Space
    ' Main raytracing routine

    Dim Ray(7)
    Dim xhHit(3), xtHit(3), xnHit(3)
    Dim xhReflect(3), xtReflect(3), xnReflect(3)
    Dim xhShadow(3), xtShadow(3), xnShadow(3)

    ' Main loop: one calculation per pixel
    For Ny% = 0 To ScrHeight - 1
        For Nx% = 0 To ScrWidth - 1

            ' Probe the objects with the pixel ray
            Probe% = 1: AntiAlias = 0
            HitAliasCheck:
            Ray.Pixel Nx%, Ny%, Ray(): ShapeHit% = 0
            Probe.Ray Ray(), ShapeHit%, Dist, xhHit(), xtHit(), xnHit(), FlagHit%

            ' AntiAlias checking: iterate until two similar hits in case of a change
            If shape(ShapeHit%).Degree = 4 Or shape(LastShape%).Degree = 4 Then
                If Probe% = 1 And (ShapeHit% <> LastShape% Or FlagHit% <> LastFlag%) Then
                    Probe% = 2: AntiAlias = 1
                    Hit1% = ShapeHit%: Flag1% = FlagHit%
                    GoTo HitAliasCheck:
                End If
                If Probe% = 2 And (ShapeHit% <> Hit1% Or FlagHit% <> Flag1%) Then
                    Probe% = 2: AntiAlias = 0
                    Hit1% = ShapeHit%: Flag1% = FlagHit%
                    GoTo HitAliasCheck:
                End If
            End If
            LastShape% = ShapeHit%: LastFlag% = FlagHit%

            ' Calculate the pixel colour in case a shape is hit
            If FlagHit% <> 0 Then

                ' Check if a shadow must be applied
                Probe% = 1: AntiAlias = 0
                d1 = xtHit(1) - xLight(1)
                d2 = xtHit(2) - xLight(2)
                d3 = xtHit(3) - xLight(3)
                DistToLight = Module(d1, d2, d3)
                ShadowAliasCheck:
                Ray.Shadow xtHit(), Ray(): ShapeShadow% = ShapeHit%
                Probe.Ray Ray(), ShapeShadow%, DistShadow, xhShadow(), xtShadow(), xnShadow(), FlagShadow%
                'IF DistShadow > DistToLight THEN FlagShadow% = 0
                If DistShadow > 30 Then FlagShadow% = 0

                ' AntiAlias checking: iterate until two similar hits in case of a change
                If shape(ShapeShadow%).Degree = 4 Or shape(LastShadow%).Degree = 4 Then
                    If Probe% = 1 And FlagShadow% <> LastShadow% Then
                        Probe% = 2: AntiAlias = 0
                        Flag1% = FlagShadow%
                        GoTo ShadowAliasCheck:
                    End If
                    If Probe% = 2 And FlagShadow% <> Flag1% Then
                        Probe% = 2: AntiAlias = 1
                        Flag1% = FlagShadow%
                        GoTo ShadowAliasCheck:
                    End If
                End If
                LastShadow% = FlagShadow%

                If FlagShadow% <> 0 Then
                    Shadow = cShadow / (1 + .5 * DistShadow) * shape(ShapeHit%).Brilliance
                Else
                    Shadow = shape(ShapeHit%).Brilliance
                End If

                ' Calculate the pixel colour without reflection
                HSVtoRGB Texture.Object!(ShapeHit%, FlagHit%, xhHit(), xnHit()), shape(ShapeHit%).SatRef, Shadow * Light.Specular(xtHit(), xnHit())

                ' Check if a reflection colour must be added
                If shape(ShapeHit%).Reflect <> 0 Then

                    Probe% = 1: AntiAlias = 0
                    ReflectAliasCheck:
                    Ray.Reflect xtHit!(), xnHit!(), Ray(): ShapeReflect% = ShapeHit%
                    Probe.Ray Ray(), ShapeReflect%, DistReflect, xhReflect(), xtReflect(), xnReflect(), FlagReflect%

                    ' AntiAlias checking: iterate until two similar hits in case of a change
                    If shape(ShapeReflect%).Degree = 4 Or shape(LastReflect%).Degree = 4 Then
                        If Probe% = 1 And FlagReflect% <> LastReflect% Then
                            Probe% = 2: AntiAlias = 0
                            Flag1% = FlagReflect%
                            GoTo ReflectAliasCheck:
                        End If
                        If Probe% = 2 And FlagReflect% <> Flag1% Then
                            Probe% = 2: AntiAlias = 1
                            Flag1% = FlagReflect%
                            GoTo ReflectAliasCheck:
                        End If
                    End If
                    LastReflect% = FlagReflect%

                    ' Merge the two colours in case of reflection
                    If FlagReflect% <> 0 Then
                        r1% = Red%: g1% = Green%: B1% = Blue%
                        HSVtoRGB Texture.Object!(ShapeReflect%, FlagReflect%, xtReflect(), xnReflect()), shape(ShapeReflect%).SatRef, shape(ShapeReflect%).Brilliance * Light.Specular(xtReflect(), xnReflect())
                        r2% = Red%: g2% = Green%: B2% = Blue%
                        Reflection = shape(ShapeHit%).Reflect / (1 + .15 * DistReflect)
                        Red% = r1% + Reflection * r2%
                        Green% = g1% + Reflection * g2%
                        Blue% = B1% + Reflection * B2%
                        If Red% > 255 Then Red% = 255
                        If Green% > 255 Then Green% = 255
                        If Blue% > 255 Then Blue% = 255
                    End If
                End If

                ' Apply Gamma and exposure corrections if wanted
                If Gamma <> 1 Then
                    RGBtoHSV Hue, Sat, Value
                    'gamma correction
                    Value = ((Value / 100) ^ (1 / Gamma)) * 100
                    'exposure control
                    Value = (1 - Exp(-Value * .01 * Aperture)) * 100
                    HSVtoRGB Hue, Sat, Value
                End If

                ' Plot the point...
                tt& = _RGB(Red%, Green%, Blue%)
                PReset (Nx%, Ny%), tt&
            End If

        Next Nx%
        _Display
        If Len(InKey$) Then Exit Sub 'agv. to be able to kill the rendering
    Next Ny%


End Sub

Sub Equa1 (A1#, a0#, x1#, Err1%, nSol1%)
    ' First degree equation resolution
    If A1# = 0 And a0# = 0 Then Err1% = 1: nSol1% = 0: Exit Sub
    If A1# = 0 And a0# <> 0 Then Err1% = 0: nSol1% = 0: Exit Sub
    If A1# <> 0 Then Err1% = 0: x1# = -a0# / A1#: nSol1% = 1: Exit Sub
End Sub

Sub Equa2 (A2#, A1#, a0#, x1#, x2#, Err2%, nSol2%)
    ' Second degree equation resolution - Discriminant method
    Err2% = 0
    If A2# = 0 And A1# = 0 And a0# = 0 Then Err2% = 1: nSol2% = 0: Exit Sub
    If A2# = 0 And A1# = 0 And a0# <> 0 Then nSol2% = 0: Exit Sub

    If A2# = 0 Then
        Equa1 A1#, a0#, x1#, Err2%, nSol2%
    Else
        Delta2# = A1# * A1# - 4 * A2# * a0#
        If Delta2# >= 0 Then
            nSol2% = 2
            x1# = -(A1# + Sqr(Delta2#)) * .5 / A2#
            x2# = (Sqr(Delta2#) - A1#) * .5 / A2#
        Else
            nSol2% = 0: Exit Sub
        End If
    End If
End Sub

Sub Equa3 (a30#, a20#, a10#, a00#, x1#, x2#, x3#, Err3%, nSol3%)
    ' Third degree equation resolution - Cardan method
    Err3% = 0

    If a30# = 0 And a20# = 0 And a10# = 0 And a00# = 0 Then
        nSol3% = 0
        Err3% = 1
        Exit Sub
    End If

    If a30# = 0 And a20# = 0 And a10# = 0 And a00# <> 0 Then
        nSol3% = 0
        Exit Sub
    End If

    If a30# = 0 And a20# = 0 Then
        Equa1 a10#, a00#, x1#, Err3%, nSol3%
        Exit Sub
    End If

    If a30# = 0 Then
        Equa2 a20#, a10#, a00#, x1#, x2#, Err3%, nSol3%
        Exit Sub
    End If

    A3# = 1
    A2# = a20# / a30#
    A1# = a10# / a30#
    a0# = a00# / a30#

    p# = A1# - 1 / 3 * (A2#) ^ 2
    q# = 2 / 27 * (A2#) ^ 3 + a0# - A2# / 3 * A1#
    Delta3# = q# ^ 2 + 4 / 27 * p# ^ 3

    If Delta3# > 0 Then
        x1# = Sqr3(-q# / 2 + Sqr(Delta3#) / 2) + Sqr3(-q# / 2 - Sqr(Delta3#) / 2) - A2# / 3
        nSol3% = 1
        Exit Sub
    End If

    nSol3% = 3

    If Delta3# = 0 Then
        If p# <> 0 Then
            x1# = q# / p# * 3 - A2# / 3
            x2# = -Sqr(-p# / 3) - A2# / 3
            x3# = x2#
            If x1# > x2# Then Swap x1#, x2#
            If x1# > x3# Then Swap x1#, x3#
            If x2# > x3# Then Swap x2#, x3#
        Else
            x1# = -A2# / 3: x2# = x1#: x3# = x1#
        End If
        Exit Sub
    End If

    If Delta3# < 0 Then
        r# = Sqr(-p# ^ 3 / 27)
        CosTeta3# = -q# / 2 / r#
        Teta3# = Acos(CosTeta3#)

        x1# = 2 * Sqr(-p# / 3) * Cos(Teta3# / 3) - A2# / 3
        x2# = 2 * Sqr(-p# / 3) * Cos(Teta3# / 3 + 2 * Pi / 3) - A2# / 3
        x3# = 2 * Sqr(-p# / 3) * Cos(Teta3# / 3 + 4 * Pi / 3) - A2# / 3
        If x1# > x2# Then Swap x1#, x2#
        If x1# > x3# Then Swap x1#, x3#
        If x2# > x3# Then Swap x2#, x3#
    End If

End Sub

Sub Equa4 (a40#, a30#, a20#, a10#, a00#, x1#, x2#, x3#, x4#, Err4%, nSol4%)
    ' 4th degree polynomial equations - Ferrari method

    ' If a40#=0, the equation is a cubic
    If a40# = 0 Then
        Equa3 a30#, a20#, a10#, a00#, x1#, x2#, x3#, Err4%, nSol4%
        Exit Sub
    End If

    ' Change the coefficients so that a4#=1
    a4# = 1
    A3# = a30# / a40#
    A2# = a20# / a40#
    A1# = a10# / a40#
    a0# = a00# / a40#

    ' The equation to be solved is now
    ' x^4 + a3*x^3 + a2*x^2 + a1*x + a0 = 0  (1)

    ' Change variable to remove the cubic term
    ' X = x + a3#/4
    ' The original equation becomes:
    ' X^4 + A*X^2 + B*X + C = 0     (2)

    a# = A2# - 3 / 8 * A3# ^ 2
    b# = A3# ^ 3 / 8 - A3# * A2# / 2 + A1#
    C# = A3# ^ 2 * A2# / 16 - 3 / 256 * A3# ^ 4 - A1# * A3# / 4 + a0#

    ' If B = 0 then (2) is a bi-quadric equation:
    If b# = 0 Then
        Equa2 1, a#, C#, x21#, x22#, Err2%, nSol2%
        If x21# < 0 And x22# < 0 Then nSol4% = 0: Exit Sub
        If x21# >= 0 And x22# >= 0 Then
            nSol4% = 4
            x1# = Sqr(x21#) - A3# / 4
            x2# = -Sqr(x21#) - A3# / 4
            x3# = Sqr(x22#) - A3# / 4
            x4# = -Sqr(x22#) - A3# / 4
            Exit Sub
        End If
        If x21# >= 0 And x22# < 0 Then
            nSol4% = 2
            x1# = Sqr(x21#) - A3# / 4
            x2# = -Sqr(x21#) - A3# / 4
            x3# = 0
            x4# = 0
            Exit Sub
        End If
        If x21# < 0 And x22# >= 0 Then
            nSol4% = 2
            x1# = Sqr(x22#) - A3# / 4
            x2# = -Sqr(x22#) - A3# / 4
            x3# = 0
            x4# = 0
            Exit Sub
        End If

    End If

    ' If B<>0, we must introduce an auxiliary variable "u".
    ' Let's calculate (X^2 + u/2)^2:
    ' (X^2 + u/2)^2 = (u-A)*X^2 - B*X + u^2/4 - C  (3)

    ' We want (3) to be equivalent to (2) for any "u"
    ' That means the second term of (3) is a double-root second degree polynom
    ' This condition is obtained if its discriminant is nul

    ' That leads to the cubic solvant equation:
    ' u^3 - A * u^2 - 4 * C * u + (4*A*C - B^2) = 0

    p# = -a#
    q# = -4 * C#
    r# = 4 * a# * C# - b# * b#
    Equa3 1, p#, q#, r#, u1#, u2#, u3#, Err3%, nSol3%

    ' This cubic equation has a root U0, with U0 > A
    ' If such a root does't exist, then B = 0, which means equation (2)
    ' is a bi-quadric equation.

    If u1# > a# Then
        u0# = u1#
    Else
        If u2# > a# Then
            u0# = u2#
        Else
            If u3# > a# Then
                u0# = u3#
            Else
                nSol4% = 0
            End If
        End If
    End If

    ' Equation (3) can now be re-written:
    ' (X^2 + u0/2)^2 = (u0 - A) * (X - B/2/(u0-A))^2
    ' This bi-quadric equation is equivalent to the two following quadrics:

    ' X^2 + u0/2 = +SQR(u0-A) * X - B/2/SQR(u0-A)
    ' X^2 + u0/2 = -SQR(u0-A) * X + B/2/SQR(u0-A)

    SF# = Sqr(Abs(u0# - a#))
    Equa2 1, -SF#, (u0# / 2 + b# / 2 / SF#), xa1#, xa2#, Err2a%, nSol2a%
    Equa2 1, SF#, (u0# / 2 - b# / 2 / SF#), xb1#, xb2#, Err2b%, nSol2b%

    If nSol2a% = 0 And nSol2b% = 0 Then nSol4% = 0: Exit Sub

    If nSol2a% = 2 And nSol2b% = 0 Then
        nSol4% = 2
        x1# = xa1# - A3# / 4
        x2# = xa2# - A3# / 4
        x3# = 0
        x4# = 0
        Exit Sub
    End If

    If nSol2a% = 0 And nSol2b% = 2 Then
        nSol4% = 2
        x1# = xb1# - A3# / 4
        x2# = xb2# - A3# / 4
        x3# = 0
        x4# = 0
        Exit Sub
    End If

    If nSol2a% = 2 And nSol2b% = 2 Then
        nSol4% = 4
        x1# = xa1# - A3# / 4
        x2# = xa2# - A3# / 4
        x3# = xb1# - A3# / 4
        x4# = xb2# - A3# / 4
        Exit Sub
    End If

End Sub

Sub Init.Cubic (Id%, Type$, HueRef, x0!, y0!, z0!, dRef1, dRef2, Alpha, Beta)
    ' The power of three... somehow dangerous!
    ' I'm not sure this routine is 100% safe

    Select Case Type$
        Case "C0"
            shape(Id%).NumType = 300
            KXYZ = 0: KX3 = 1: KY3 = 1: KZ3 = 1
            KX2Y = 0: KX2Z = 0: KY2Z = 0
            KY2X = 0: KZ2X = 0: KZ2Y = 0
            KX2 = 0: KY2 = 0: KZ2 = 0
            KXY = 0: KXZ = 0: KYZ = 0
            kx = 0: ky = 0: KZ = 0
            K0 = -dRef1 ^ 3

        Case "C1"
            shape(Id%).NumType = 301
            KXYZ = 0: KX3 = 0: KY3 = 0: KZ3 = 1
            KX2Y = 0: KX2Z = 1: KY2Z = 1
            KY2X = 0: KZ2X = 0: KZ2Y = 0
            KX2 = 0: KY2 = 0: KZ2 = 0
            KXY = 0: KXZ = 0: KYZ = 0
            kx = 0: ky = 0: KZ = -dRef1 ^ 2
            K0 = 0
        Case Else: Exit Sub
    End Select

    For k% = 1 To 15: Object(Id%, k%) = 0: Next k%
    Object(Id%, 16) = KXYZ
    Object(Id%, 17) = KX3: Object(Id%, 18) = KY3: Object(Id%, 19) = KZ3
    Object(Id%, 20) = KX2Y: Object(Id%, 21) = KX2Z: Object(Id%, 22) = KY2Z
    Object(Id%, 23) = KY2X: Object(Id%, 24) = KZ2X: Object(Id%, 25) = KZ2Y
    Object(Id%, 26) = KX2: Object(Id%, 27) = KY2: Object(Id%, 28) = KZ2
    Object(Id%, 29) = KXY: Object(Id%, 30) = KYZ: Object(Id%, 31) = KXZ
    Object(Id%, 32) = kx: Object(Id%, 33) = ky: Object(Id%, 34) = KZ
    Object(Id%, 35) = K0

    shape(Id%).Degree = 3
    shape(Id%).ShapeType = Type$
    shape(Id%).HueRef = HueRef: shape(Id%).SatRef = 200 / 3
    shape(Id%).dRef1 = dRef1: shape(Id%).dRef2 = dRef2
    shape(Id%).x0 = x0!: shape(Id%).y0 = y0!: shape(Id%).z0 = z0!
    shape(Id%).Brilliance = 1
    shape(Id%).Reflect = cReflect

    Matrix(Id%, 1, 1) = CosD(Beta) * CosD(Alpha)
    Matrix(Id%, 1, 2) = CosD(Beta) * SinD(Alpha)
    Matrix(Id%, 1, 3) = SinD(Beta)
    Matrix(Id%, 2, 1) = -SinD(Alpha)
    Matrix(Id%, 2, 2) = CosD(Alpha)
    Matrix(Id%, 2, 3) = 0
    Matrix(Id%, 3, 1) = -SinD(Beta) * CosD(Alpha)
    Matrix(Id%, 3, 2) = -SinD(Beta) * SinD(Alpha)
    Matrix(Id%, 3, 3) = CosD(Beta)

End Sub

Sub Init.Example1
    ' Geometry
    Init.Space 0, 15, 5, .8, 30, 60
    Init.Object "Fl", 240, 0, 0, 0, .15, .15, 0, 0
    shape(1).Reflect = .1
    Init.Object "H1", 200, 1, 1.2, 1, .8, .8, 20, -50
    Init.Object "Sp", 120, -1, -2.3, .5, 1.8, 1.6, -10, -20
    Init.Object "Q4", 30, 3.5, 0, .5, .4, .2, 30, -20
    shape(1).Texture = 17
End Sub

Sub Init.Example10
    'nephroid cup
    Init.Space 15, 35, 2.5, .8, 55, 30
    'Init.Space 0, 90, 2.5, .8, 55, 30

    ' Table surface
    Init.Object "Fl", 300, 0, 0, 0, 2, 2, 0, 0
    shape(1).Reflect = .9
    shape(1).Texture = 12

    ' External surface of the cup
    Init.Object "Cp", 60, .3, 0, .9, .8, 1.45, 30, 0
    shape(2).Reflect = .1
    shape(2).Texture = 11

    ' Internal surface of the cup
    Init.Object "Cp", 150, .3, 0, .9, .75, 1.45, 0, 0
    shape(3).Reflect = 0

    ' Close the cup with a torus ring
    Init.Object "T0", 240, .3, 0, .9, (.8 + .75) / 2, .025, 0, 0
    shape(4).Reflect = 0

    ' Coffee surface
    'x^2+y^2=-(dRef1/dRef2)^2*z^2 + dRef1^2
    dRef1 = .75: dRef2 = 1.45: dZ = .1
    Radius = Sqr(dRef1 ^ 2 - (dRef1 / dRef2) ^ 2 * dZ ^ 2)
    Init.Object "D1", 30, .3, 0, .8, Radius, Radius, 0, 0
    shape(5).Reflect = 0
    shape(5).SatRef = 50
    shape(5).Texture = 32

    ' Torus for the cup handle
    Init.Object "T0", 240, .3, .7, .4, .3, .05, 0, 90
    shape(6).Reflect = 0
End Sub

Sub Init.Example2
    ' Birthday
    Init.Space 0, 30, 5, .8, 30, 60
    Init.Object "Fl", 240, 0, 0, 0, .15, .15, 0, 0
    shape(1).Texture = 17

    r! = 2
    For k% = 1 To 16
        a! = k% / 16 * 2 * Pi
        Init.Object "Bl", k% * 24, r! * Cos(a!), r! * Sin(a!), 1.5, .3, .3, 0, 0
    Next k%
    r! = 2.5
    For k% = 1 To 16
        a! = k% / 16 * 2 * Pi
        Init.Object "Cy", 180 + k% * 24, r! * Cos(a!), r! * Sin(a!), .4, .3, .8, 0, 0
    Next k%


End Sub

Sub Init.Example3
    ' ADN Helix
    Init.Space 0, 40, 5, .8, 30, 60
    Dim u(3), v(3), W(3)

    Init.Object "Fl", 90, 0, 0, 0, .3, .3, 0, 0
    shape(1).Texture = 17
    r! = 1.6
    Alpha = Pi / 4: Beta = Pi / 3
    u(1) = Cos(Alpha) * Cos(Beta): v(1) = -Sin(Alpha): W(1) = -Sin(Beta) * Cos(Alpha)
    u(2) = Sin(Alpha) * Cos(Beta): v(2) = Cos(Alpha): W(2) = -Sin(Beta) * Sin(Alpha)
    u(3) = Sin(Beta): v(3) = 0: W(3) = Cos(Beta)

    For k% = 1 To 48
        a! = k% / 16 * 2 * Pi
        x0r = r! * Cos(a!)
        y0r = r! * Sin(a!)
        z0r = -4 + k% * .15
        x0! = u(1) * x0r + v(1) * y0r + W(1) * z0r
        y0! = u(2) * x0r + v(2) * y0r + W(2) * z0r
        z0! = u(3) * x0r + v(3) * y0r + W(3) * z0r
        z0! = z0! + r! / 3
        Init.Object "Bl", k% * 24, x0!, y0!, z0!, .4, .4, 0, 0
    Next k%

End Sub

Sub Init.Example4
    ' Perlin Ring
    Init.Space 23, 30, 5, .8, 30, 60
    cShadow = .7: cReflect = .2

    Init.Object "T0", 120, 0, 0, .8, 2, .4, 0, 0
    shape(1).Texture = 16

    For k% = 1 To 4
        a! = k% / 2 * Pi
        Init.Object "Bl", k% * 90, 2 * Cos(a!), 2 * Sin(a!), .8, .8, .8, 0, 0
    Next k%
    shape(2).Texture = 13
    shape(3).Texture = 14
    shape(4).Texture = 15
    shape(5).Texture = 12

    Init.Object "Fl", 0, 0, 0, 0, 1, 1, 0, 0
    shape(6).Texture = 13
End Sub

Sub Init.Example5
    ' Stormy Sea
    Init.Space 5, 10, 5, .6, 40, 30
    Init.Object "Sk", 260, 0, 0, 20, 0, 0, 0, 0
    Init.Object "Fl", 190, 0, 0, 0, 1, 1, 0, 0
    Init.Object "HC", 0, 0, 0, .7, .5, .2, -30, -15
    shape(2).Texture = 30
    shape(3).Texture = 20
End Sub

Sub Init.Example6
    ' Dolmen
    Init.Space 10, 8, 4.9, .6, 40, 30
    Init.Object "Fl", 120, 0, 0, 0, 2, 2, 0, 0
    shape(1).Texture = 30
    shape(1).Reflect = .3
    Init.Object "Sk", 190, 0, 0, 20, 0, 0, 0, 0

    Scale = 1.25
    Init.Object "T1", 120, 0, 0, 1.6 * Scale, 1.2 * Scale, .2 * Scale, 0, 0
    Init.Object "El", 120, Scale, Scale, .5 * Scale, .3 * Scale, .7 * Scale, 0, 0
    Init.Object "El", 120, Scale, -Scale, .5 * Scale, .3 * Scale, .7 * Scale, 0, 0
    Init.Object "El", 120, -Scale, Scale, .5 * Scale, .3 * Scale, .7 * Scale, 0, 0
    Init.Object "El", 120, -Scale, -Scale, .5 * Scale, .3 * Scale, .7 * Scale, 0, 0

    Scale = 2
    Init.Object "El", 120, 2.5, -1.4, .5 * Scale, .2 * Scale, .7 * Scale, 0, 0
    Scale = 1.25
    Init.Object "El", 120, 3, 2, .5 * Scale, .3 * Scale, .7 * Scale, 0, 0

    For k% = 3 To NbObjects%
        shape(k%).SatRef = .5
        shape(k%).Bump = 1
        shape(k%).Reflect = .2
        shape(k%).Brilliance = .4
    Next k%

End Sub

Sub Init.Example7
    ' MandelChrist
    ' Some comments are necessary
    ' This pic does not want to hurt nore shock anybody
    ' The one who would laugh or even smile after it can only be a stupid *#\
    ' Let's say I believe in Jesus, cause he was a man of faith in other men
    ' I just don't want this world to end in a sea of blood
    ' Hope this pic can help
    ' .

    Init.Space -30, 20, 5, .8, 25, 40
    Init.Object "Gm", 300, .4, 0, .81, 0, 0, 0, 0
    shape(1).Reflect = 0
    shape(1).Brilliance = .5
    Init.Object "Sk", 240, 0, 0, 20, 0, 0, 0, 0
    Init.Object "Fl", 0, 0, 0, 0, 1, 1, 0, 0
    shape(3).Texture = 30
    CrossHue = 60
    Init.Object "P2", CrossHue, .4, 0, 0, .4, 2.5, 0, 0 'Vertical bar, front
    Init.Object "P2", CrossHue, .4, 0, 1, 2, .4, 0, 0 'Horizontal bar, front
    Init.Object "P2", CrossHue, -.4, 0, 0, .4, 2.5, 0, 0 'Vertical bar, rear
    Init.Object "P2", CrossHue, -.4, 0, 1, 2, .4, 0, 0 'Horizontal bar, rear
    Init.Object "P1", CrossHue, 0, 0, 1.4, .4, 2, 0, 0 'Horizontal bar, top
    Init.Object "P1", CrossHue, 0, 0, .6, .4, 2, 0, 0 'Horizontal bar, bottom
    Init.Object "P3", CrossHue, 0, -.4, 0, .4, 2.5, 0, 0 'Vertical bar, left
    Init.Object "P3", CrossHue, 0, .4, 0, .4, 2.5, 0, 0 'Vertical bar, right
    Init.Object "P3", CrossHue, 0, -2, 1, .4, .4, 0, 0 'Horizontal bar, left
    Init.Object "P3", CrossHue, 0, 2, 1, .4, .4, 0, 0 'Horizontal bar, right

    For k% = 4 To NbObjects%
        shape(k%).Bump = 2
        shape(k%).Reflect = 0
        shape(k%).Brilliance = .4
        shape(k%).SatRef = 200 / 3
    Next k%
End Sub

Sub Init.Example8
    ' Bike Mirrors
    Init.Space 0, 15, 5, .8, 40, 30
    Init.Object "C0", 30, 0, 0, 0, 1, 8, -20, -20
    shape(1).Texture = 12
    Init.Object "Q4", 240, 1.6, 0, .8, 1, .5, 20, -30
    shape(2).Texture = 16
End Sub

Sub Init.Example9
    ' Water Effect
    Init.Space 0, 15, 5, .8, 30, 60
    Init.Object "Fl", 240, 0, 0, 0, 1, 1, 0, 0
    Init.Object "Bl", 240, 0, 0, 1, 1, 1, 0, 0
    shape(1).Bump = 3
End Sub

Sub Init.Object (Type$, HueRef, Px0!, Py0!, Pz0!, dRef1, dRef2, Alpha, Beta)
    ' Select the proper routine to create the object
    ' Increment number of objects at each call

    NbObjects% = NbObjects% + 1
    Select Case Type$
        Case "Fl", "P1", "P2", "P3", "D1", "D2", "D3", "Sk", "Se", "Gm"
            Init.Plan NbObjects%, Type$, HueRef, Px0!, Py0!, Pz0!, dRef1, dRef2, Alpha, Beta
        Case "Bl", "Sp", "El", "Cy", "Cn", "Co", "H1", "H2", "Hp", "Pa", "Cp"
            Init.Quadric NbObjects%, Type$, HueRef, Px0!, Py0!, Pz0!, dRef1, dRef2, Alpha, Beta
        Case "C0", "C1"
            Init.Cubic NbObjects%, Type$, HueRef, Px0!, Py0!, Pz0!, dRef1, dRef2, Alpha, Beta
        Case "T0", "T1", "HC", "Q0", "Q1", "Q2", "Q3", "Q4", "S2", "S3", "L0"
            Init.Quartic NbObjects%, Type$, HueRef, Px0!, Py0!, Pz0!, dRef1, dRef2, Alpha, Beta

        Case Else

            Print "Non valid object type: "; Type$
            Sleep
            EndProg
    End Select
End Sub

Sub Init.Perlin
    ' Inits the Noise% array for 2D perlin noises
    ' The pattern is perfectly tileable
    ' Based on squares division technique
    ' Fast, but can probably be optimised more:
    ' 2^n+1 can probably be replaced by plain 2^n via AND or MOD tricks...

    Randomize Timer: Seed& = 103476 * Rnd ' This number fell from my keyboard :-)
    Persistence! = 1.5: Span% = 2 ^ octaves% + 1

    For k% = 0 To octaves% + 1
        Amplitude!(k%) = (1 / Persistence!) ^ k% ' You can try different logics!
    Next k%

    ' Remark: the Init corners and edges lines make the pattern tileable

    ' Init the corners
    Noise%(1, 1) = 128
    Noise%(1, Span%) = 128
    Noise%(Span%, 1) = 128
    Noise%(Span%, Span%) = 128

    ' Init the edges
    For Rank% = 1 To octaves%

        Grid% = 2 ^ (octaves% - Rank% + 1)
        nStep% = 2 ^ (Rank% - 1) + 1

        For kx% = 1 To nStep% - 1
            x% = (kx% - 1) * Grid% + 1: y% = 1
            Alt% = (Noise%(x%, y%) + Noise%(x% + Grid%, y%)) / 2
            zNew% = Int(Alt% * (1 + (Rnd - .5) * Amplitude!(Rank%)))
            Noise%(x% + Grid% / 2, 1) = zNew%
            Noise%(x% + Grid% / 2, Span%) = zNew%
        Next kx%

        For ky% = 1 To nStep% - 1
            x% = 1: y% = (ky% - 1) * Grid% + 1
            Alt% = (Noise%(x%, y%) + Noise%(x%, y% + Grid%)) / 2
            zNew% = Int(Alt% * (1 + (Rnd - .5) * Amplitude!(Rank%)))
            Noise%(1, y% + Grid% / 2) = zNew%
            Noise%(Span%, y% + Grid% / 2) = zNew%
        Next ky%

    Next Rank%

    ' Now, the routine really becomes random...
    ' Fill the array
    For Rank% = 1 To octaves%

        Grid% = 2 ^ (octaves% - Rank% + 1)
        nStep% = 2 ^ (Rank% - 1) + 1

        For kx% = 1 To nStep% - 1
            x% = (kx% - 1) * Grid% + 1
            For ky% = 1 To nStep% - 1
                y% = (ky% - 1) * Grid% + 1

                Alt% = (Noise%(x%, y%) + Noise%(x% + Grid%, y%) + Noise%(x%, y% + Grid%) + Noise%(x% + Grid%, y% + Grid%)) / 4
                Noise%(x% + Grid% / 2, y% + Grid% / 2) = Int(Alt% * (1 + (Rnd - .5) * Amplitude!(Rank%)))

                Alt% = (Noise%(x%, y%) + Noise%(x% + Grid%, y%)) / 2
                If y% <> 1 Then Noise%(x% + Grid% / 2, y%) = Int(Alt% * (1 + (Rnd - .5) * Amplitude!(Rank%)))

                Alt% = (Noise%(x%, y%) + Noise%(x%, y% + Grid%)) / 2
                If x% <> 1 Then Noise%(x%, y% + Grid% / 2) = Int(Alt% * (1 + (Rnd - .5) * Amplitude!(Rank%)))

                Alt% = (Noise%(x% + Grid%, y%) + Noise%(x% + Grid%, y% + Grid%)) / 2
                If (x% + Grid%) <> Span% Then Noise%(x% + Grid%, y% + Grid% / 2) = Int(Alt% * (1 + (Rnd - .5) * Amplitude!(Rank%)))

                Alt% = (Noise%(x%, y% + Grid%) + Noise%(x% + Grid%, y% + Grid%)) / 2
                If (y% + Grid%) <> Span% Then Noise%(x% + Grid% / 2, y% + Grid%) = Int(Alt% * (1 + (Rnd - .5) * Amplitude!(Rank%)))

            Next ky%
        Next kx%

    Next Rank%


End Sub

Sub Init.Plan (Id%, Type$, HueRef, Px0!, Py0!, Pz0!, dRef1, dRef2, Alpha, Beta)
    ' A plan is nothing but 4 figures, i.e. a normal vector and a shift
    ' I found it handy to have vertical and horizontal plans,
    ' and then to rotate them at will
    ' Check an application at Example 8...
    ' and please don't be shocked by the result !

    Select Case Type$
        Case "Fl": kx = 0: ky = 0: KZ = 1: K0 = 0: shape(Id%).NumType = 100
        Case "P1": kx = 0: ky = 0: KZ = 1: K0 = 0: shape(Id%).NumType = 101
        Case "P2": kx = 1: ky = 0: KZ = 0: K0 = 0: shape(Id%).NumType = 102
        Case "P3": kx = 0: ky = 1: KZ = 0: K0 = 0: shape(Id%).NumType = 103
        Case "D1": kx = 0: ky = 0: KZ = 1: K0 = 0: shape(Id%).NumType = 104
        Case "D2": kx = 1: ky = 0: KZ = 0: K0 = 0: shape(Id%).NumType = 105
        Case "D3": kx = 0: ky = 1: KZ = 0: K0 = 0: shape(Id%).NumType = 106
        Case "Sk": kx = 0: ky = 0: KZ = 1: K0 = 0: shape(Id%).NumType = 110
        Case "Gm": kx = 1: ky = 0: KZ = 0: K0 = 0: shape(Id%).NumType = 120

        Case Else
            Exit Sub
    End Select

    For k% = 1 To 31: Object(Id%, k%) = 0: Next k%
    Object(Id%, 32) = kx: Object(Id%, 33) = ky: Object(Id%, 34) = KZ: Object(Id%, 35) = K0

    shape(Id%).Degree = 1
    shape(Id%).ShapeType = Type$
    shape(Id%).HueRef = HueRef: shape(Id%).SatRef = 200 / 3
    shape(Id%).dRef1 = dRef1: shape(Id%).dRef2 = dRef2
    shape(Id%).x0 = Px0!: shape(Id%).y0 = Py0!: shape(Id%).z0 = Pz0!
    shape(Id%).Brilliance = 1
    shape(Id%).Reflect = cReflect

    ' Sky is such a strange buddy...
    If shape(Id%).NumType = 110 Then shape(Id%).Reflect = 0
    If shape(Id%).NumType = 110 Then shape(Id%).Texture = 31

    ' Now you can rotate the plan:
    Matrix(Id%, 1, 1) = CosD(Beta) * CosD(Alpha)
    Matrix(Id%, 1, 2) = CosD(Beta) * SinD(Alpha)
    Matrix(Id%, 1, 3) = SinD(Beta)
    Matrix(Id%, 2, 1) = -SinD(Alpha)
    Matrix(Id%, 2, 2) = CosD(Alpha)
    Matrix(Id%, 2, 3) = 0
    Matrix(Id%, 3, 1) = -SinD(Beta) * CosD(Alpha)
    Matrix(Id%, 3, 2) = -SinD(Beta) * SinD(Alpha)
    Matrix(Id%, 3, 3) = CosD(Beta)
End Sub

Sub Init.Quadric (Id%, Type$, HueRef, x0!, y0!, z0!, dRef1, dRef2, Alpha, Beta)
    ' Well, that's spheres, ellipsoids, hyperboloids, and all this kind of things!

    Select Case Type$

        Case "Bl" ' Plain Ball
            shape(Id%).NumType = 200
            KX2 = 1: KY2 = 1: KZ2 = 1
            KXY = 0: KYZ = 0: KXZ = 0
            kx = 0: ky = 0: KZ = 0
            K0 = -dRef1 ^ 2

        Case "Sp" ' Sphere
            shape(Id%).NumType = 201
            KX2 = 1: KY2 = 1: KZ2 = 1
            KXY = 0: KYZ = 0: KXZ = 0
            kx = 0: ky = 0: KZ = 0
            K0 = -dRef1 ^ 2

        Case "El" ' Ellipso‹d
            shape(Id%).NumType = 202
            KX2 = 1: KY2 = 1: KZ2 = (dRef1 / dRef2) ^ 2
            KXY = 0: KYZ = 0: KXZ = 0
            kx = 0: ky = 0: KZ = 0
            K0 = -dRef1 ^ 2

        Case "Cp" ' Cup
            shape(Id%).NumType = 210
            KX2 = 1: KY2 = 1: KZ2 = (dRef1 / dRef2) ^ 2
            KXY = 0: KYZ = 0: KXZ = 0
            kx = 0: ky = 0: KZ = 0
            K0 = -dRef1 ^ 2

        Case "Cy" ' Cylinder
            shape(Id%).NumType = 203
            KX2 = 1: KY2 = 1: KZ2 = 0
            KXY = 0: KYZ = 0: KXZ = 0
            kx = 0: ky = 0: KZ = 0
            K0 = -dRef1 ^ 2

        Case "Co" ' Cone
            shape(Id%).NumType = 204
            KX2 = 1: KY2 = 1: KZ2 = -(dRef1 / dRef2) ^ 2
            KXY = 0: KYZ = 0: KXZ = 0
            kx = 0: ky = 0: KZ = 0
            K0 = 0

        Case "Cn" ' Half-Cone
            shape(Id%).NumType = 205
            KX2 = 1: KY2 = 1: KZ2 = -(dRef1 / dRef2) ^ 2
            KXY = 0: KYZ = 0: KXZ = 0
            kx = 0: ky = 0: KZ = 0
            K0 = 0

        Case "Hp" ' Hyperboloid paraboloid
            shape(Id%).NumType = 206
            KX2 = 1: KY2 = -1: KZ2 = 0
            KXY = 0: KYZ = 0: KXZ = 0
            kx = 0: ky = 0: KZ = 2 * dRef1
            K0 = 0

        Case "H1" ' Hyperbolo‹d, one sheet
            shape(Id%).NumType = 207
            KX2 = 1: KY2 = 1: KZ2 = -(dRef1 / dRef2) ^ 2
            KXY = 0: KYZ = 0: KXZ = 0
            kx = 0: ky = 0: KZ = 0
            K0 = -dRef1 ^ 2

        Case "H2" ' Hyperbolo‹d, two sheets
            shape(Id%).NumType = 208
            KX2 = 1: KY2 = 1: KZ2 = -(dRef1 / dRef2) ^ 2
            KXY = 0: KYZ = 0: KXZ = 0
            kx = 0: ky = 0: KZ = 0
            K0 = dRef1 ^ 2 / 4

        Case "Pa" ' Parabolo‹d
            shape(Id%).NumType = 209
            KX2 = 1: KY2 = 1: KZ2 = 0
            KXY = 0: KYZ = 0: KXZ = 0
            kx = 0: ky = 0: KZ = -dRef2
            K0 = dRef1

        Case Else
            Exit Sub
    End Select

    For k% = 1 To 25: Object(Id%, k%) = 0: Next k%
    Object(Id%, 26) = KX2: Object(Id%, 27) = KY2: Object(Id%, 28) = KZ2
    Object(Id%, 29) = KXY: Object(Id%, 30) = KYZ: Object(Id%, 31) = KXZ
    Object(Id%, 32) = kx: Object(Id%, 33) = ky: Object(Id%, 34) = KZ
    Object(Id%, 35) = K0
    shape(Id%).Degree = 2
    shape(Id%).ShapeType = Type$
    shape(Id%).HueRef = HueRef: shape(Id%).SatRef = 200 / 3
    shape(Id%).dRef1 = dRef1: shape(Id%).dRef2 = dRef2
    shape(Id%).x0 = x0!: shape(Id%).y0 = y0!: shape(Id%).z0 = z0!
    shape(Id%).Brilliance = 1
    shape(Id%).Reflect = cReflect

    Matrix(Id%, 1, 1) = CosD(Beta) * CosD(Alpha)
    Matrix(Id%, 1, 2) = CosD(Beta) * SinD(Alpha)
    Matrix(Id%, 1, 3) = SinD(Beta)
    Matrix(Id%, 2, 1) = -SinD(Alpha)
    Matrix(Id%, 2, 2) = CosD(Alpha)
    Matrix(Id%, 2, 3) = 0
    Matrix(Id%, 3, 1) = -SinD(Beta) * CosD(Alpha)
    Matrix(Id%, 3, 2) = -SinD(Beta) * SinD(Alpha)
    Matrix(Id%, 3, 3) = CosD(Beta)

End Sub

Sub Init.Quartic (Id%, Type$, HueRef, x0!, y0!, z0!, dRef1, dRef2, Alpha, Beta)
    ' That's the land of strange things...

    Select Case Type$

        Case "HC" ' Hollow Cube
            shape(Id%).NumType = 444
            KX4 = 1: KY4 = 1: KZ4 = 1
            KX2 = -10 * dRef1 ^ 2
            KY2 = KX2: KZ2 = KX2
            K0 = 40 * dRef1 ^ 4

        Case "T0" ' Classical Torus
            shape(Id%).NumType = 400
            KX4 = 1: KY4 = 1: KZ4 = 1
            KX2Y2 = 2: KX2Z2 = 2: KY2Z2 = 2
            KX2 = -2 * (dRef2 ^ 2 + dRef1 ^ 2)
            KY2 = -2 * (dRef2 ^ 2 + dRef1 ^ 2)
            KZ2 = 2 * (dRef1 ^ 2 - dRef2 ^ 2)
            K0 = (dRef1 ^ 2 - dRef2 ^ 2) ^ 2

        Case "T1" ' Bulbeous donut
            shape(Id%).NumType = 401
            KX4 = 1: KY4 = 1: KZ4 = 1
            KX2Y2 = 1: KX2Z2 = 1: KY2Z2 = 1
            KX2 = -2 * ((dRef1 / 5) ^ 2 + dRef1 ^ 2)
            KY2 = -2 * ((dRef1 / 5) ^ 2 + dRef1 ^ 2)
            KZ2 = 2 * (dRef1 ^ 2 - dRef2 ^ 2)
            K0 = (dRef1 ^ 2 - dRef2 ^ 2) ^ 2

        Case "Q0"
            shape(Id%).NumType = 410
            F1 = CQ1 * dRef2 ^ 2
            F2 = CQ2 * dRef2 / dRef1
            F3 = CQ3 * 4 * (dRef2 / dRef1) ^ 2

            KX4 = 1 + F3 ^ 2
            KY4 = 1 + F3 ^ 2
            KZ4 = 1
            KX3Y = -2 * F2 + 2 * F2 * F3
            KY3X = -2 * F2 - 2 * F2 * F3
            KZ2XY = -2 * F2
            KX2Y2 = 2 + F2 ^ 2 - 2 * F3 ^ 2
            KX2Z2 = 2
            KY2Z2 = 2
            KX2 = -2 * dRef1 ^ 2 - 2 * F1 + 2 * F1 * F3
            KY2 = -2 * dRef1 ^ 2 - 2 * F1 - 2 * F1 * F3
            KZ2 = 2 * dRef1 ^ 2 - 2 * F1
            KXY = -2 * F2 * dRef1 ^ 2 + 2 * F1 * F2
            K0 = dRef1 ^ 4 - 2 * F1 * dRef1 ^ 2 + F1 ^ 2


        Case "Q1" ' Orange slices
            shape(Id%).NumType = 411
            F1 = 0
            F2 = dRef2 / dRef1
            F3 = 4 * (dRef2 / dRef1) ^ 2

            KX4 = 1 + F3 ^ 2
            KY4 = 1 + F3 ^ 2
            KZ4 = 1
            KX3Y = -2 * F2 + 2 * F2 * F3
            KY3X = -2 * F2 - 2 * F2 * F3
            'KZ2XY = -2 * F2
            KX2Y2 = 2 + F2 ^ 2 - 2 * F3 ^ 2
            KX2Z2 = 2
            KY2Z2 = 2
            KX2 = -2 * dRef1 ^ 2 - 2 * F1 + 2 * F1 * F3
            KY2 = -2 * dRef1 ^ 2 - 2 * F1 - 2 * F1 * F3
            KZ2 = 2 * dRef1 ^ 2 - 2 * F1
            KXY = -2 * F2 * dRef1 ^ 2 + 2 * F1 * F2
            K0 = dRef1 ^ 4 - 2 * F1 * dRef1 ^ 2 + F1 ^ 2

        Case "Q2" ' Two parts cushion
            shape(Id%).NumType = 412
            F1 = dRef2 ^ 2
            F2 = dRef2 / dRef1
            F3 = 0

            KX4 = 1 + F3 ^ 2
            KY4 = 1 + F3 ^ 2
            KZ4 = 1
            KX3Y = -2 * F2 + 2 * F2 * F3
            KY3X = -2 * F2 - 2 * F2 * F3
            'KZ2XY = -2 * F2
            KX2Y2 = 2 + F2 ^ 2 - 2 * F3 ^ 2
            KX2Z2 = 2
            KY2Z2 = 2
            KX2 = -2 * dRef1 ^ 2 - 2 * F1 + 2 * F1 * F3
            KY2 = -2 * dRef1 ^ 2 - 2 * F1 - 2 * F1 * F3
            KZ2 = 2 * dRef1 ^ 2 - 2 * F1
            KXY = -2 * F2 * dRef1 ^ 2 + 2 * F1 * F2
            K0 = dRef1 ^ 4 - 2 * F1 * dRef1 ^ 2 + F1 ^ 2

        Case "Q3" ' Feminine breast
            shape(Id%).NumType = 413
            F1 = dRef2 ^ 2
            F2 = 0
            F3 = 4 * (dRef2 / dRef1) ^ 2

            KX4 = 1 + F3 ^ 2
            KY4 = 1 + F3 ^ 2
            KZ4 = 1
            KX3Y = -2 * F2 + 2 * F2 * F3
            KY3X = -2 * F2 - 2 * F2 * F3
            'KZ2XY = -2 * F2
            KX2Y2 = 2 + F2 ^ 2 - 2 * F3 ^ 2
            KX2Z2 = 2
            KY2Z2 = 2
            KX2 = -2 * dRef1 ^ 2 - 2 * F1 + 2 * F1 * F3
            KY2 = -2 * dRef1 ^ 2 - 2 * F1 - 2 * F1 * F3
            KZ2 = 2 * dRef1 ^ 2 - 2 * F1
            KXY = -2 * F2 * dRef1 ^ 2 + 2 * F1 * F2
            K0 = dRef1 ^ 4 - 2 * F1 * dRef1 ^ 2 + F1 ^ 2

        Case "Q4" ' Bike mirrors
            shape(Id%).NumType = 414
            F1 = dRef2 ^ 2
            F2 = dRef2 / dRef1
            F3 = 4 * (dRef2 / dRef1) ^ 2

            KX4 = 1 + F3 ^ 2
            KY4 = 1 + F3 ^ 2
            KZ4 = 1
            KX3Y = -2 * F2 + 2 * F2 * F3
            KY3X = -2 * F2 - 2 * F2 * F3
            'KZ2XY = -2 * F2
            KX2Y2 = 2 + F2 ^ 2 - 2 * F3 ^ 2
            KX2Z2 = 2
            KY2Z2 = 2
            KX2 = -2 * dRef1 ^ 2 - 2 * F1 + 2 * F1 * F3
            KY2 = -2 * dRef1 ^ 2 - 2 * F1 - 2 * F1 * F3
            KZ2 = 2 * dRef1 ^ 2 - 2 * F1
            KXY = -2 * F2 * dRef1 ^ 2 + 2 * F1 * F2
            K0 = dRef1 ^ 4 - 2 * F1 * dRef1 ^ 2 + F1 ^ 2


        Case "S2" ' Double sphere
            shape(Id%).NumType = 420
            KX4 = 1: KY4 = 1: KZ4 = 1
            KX2Y2 = 2: KX2Z2 = 2: KY2Z2 = 2
            KX2 = -2 * (dRef1 ^ 2 + dRef2 ^ 2)
            KY2 = 2 * (-dRef1 ^ 2 + dRef2 ^ 2)
            KZ2 = 2 * (-dRef1 ^ 2 + dRef2 ^ 2)
            K0 = (dRef1 ^ 2 - dRef2 ^ 2) ^ 2

        Case "S3" ' Pelote de laine
            shape(Id%).NumType = 421
            KX4 = 1: KY4 = 1: KZ4 = 1
            KX2Y2 = 1: KX2Z2 = 1: KY2Z2 = 1
            KX2 = -2 * (dRef1 ^ 2 + dRef2 ^ 2)
            KY2 = 2 * (-dRef1 ^ 2 + dRef2 ^ 2)
            KZ2 = 2 * (-dRef1 ^ 2 + dRef2 ^ 2)
            K0 = (dRef1 ^ 2 - dRef2 ^ 2) ^ 2

        Case "L0" ' Lemniscate
            shape(Id%).NumType = 422
            KX4 = 1: KY4 = 1: KZ4 = 1
            KX2Y2 = 2: KX2Z2 = 2: KY2Z2 = 2
            KX2 = -2 * (dRef1 ^ 2 + dRef2 ^ 2)
            KY2 = 4 * (dRef1 ^ 2 - dRef2 ^ 2)
            KZ2 = 4 * (dRef1 ^ 2 - dRef2 ^ 2)
            K0 = 0

        Case Else
            Exit Sub
    End Select


    Object(Id%, 1) = KX4: Object(Id%, 2) = KY4: Object(Id%, 3) = KZ4
    Object(Id%, 4) = KX3Y: Object(Id%, 5) = KX3Z: Object(Id%, 6) = KY3Z
    Object(Id%, 7) = KY3X: Object(Id%, 8) = KZ3X: Object(Id%, 9) = KZ3Y
    Object(Id%, 10) = KX2Y2: Object(Id%, 11) = KX2Z2: Object(Id%, 12) = KY2Z2
    Object(Id%, 13) = KX2YZ: Object(Id%, 14) = KY2XZ: Object(Id%, 15) = KZ2XY
    Object(Id%, 16) = KXYZ
    Object(Id%, 17) = KX3: Object(Id%, 18) = KY3: Object(Id%, 19) = KZ3
    Object(Id%, 20) = KX2Y: Object(Id%, 21) = KX2Z: Object(Id%, 22) = KY2Z
    Object(Id%, 23) = KY2X: Object(Id%, 24) = KZ2X: Object(Id%, 25) = KZ2Y
    Object(Id%, 26) = KX2: Object(Id%, 27) = KY2: Object(Id%, 28) = KZ2
    Object(Id%, 29) = KXY: Object(Id%, 30) = KYZ: Object(Id%, 31) = KXZ
    Object(Id%, 32) = kx: Object(Id%, 33) = ky: Object(Id%, 34) = KZ
    Object(Id%, 35) = K0

    shape(Id%).Degree = 4
    shape(Id%).ShapeType = Type$
    shape(Id%).HueRef = HueRef: shape(Id%).SatRef = 200 / 3
    shape(Id%).dRef1 = dRef1: shape(Id%).dRef2 = dRef2
    shape(Id%).x0 = x0!: shape(Id%).y0 = y0!: shape(Id%).z0 = z0!
    shape(Id%).Brilliance = 1
    shape(Id%).Reflect = cReflect

    Matrix(Id%, 1, 1) = CosD(Beta) * CosD(Alpha)
    Matrix(Id%, 1, 2) = CosD(Beta) * SinD(Alpha)
    Matrix(Id%, 1, 3) = SinD(Beta)
    Matrix(Id%, 2, 1) = -SinD(Alpha)
    Matrix(Id%, 2, 2) = CosD(Alpha)
    Matrix(Id%, 2, 3) = 0
    Matrix(Id%, 3, 1) = -SinD(Beta) * CosD(Alpha)
    Matrix(Id%, 3, 2) = -SinD(Beta) * SinD(Alpha)
    Matrix(Id%, 3, 3) = CosD(Beta)

End Sub

Sub Init.Space (AlphaCam, BetaCam, DistCam, DistScreen, aLight, bLight)

    ' Init the 2D Perlin noise array
    Init.Perlin

    ' Shadow and Reflection control coefficients
    cShadow = .3: cReflect = .5

    ' Point of light
    rLight = 10
    xLight(1) = rLight * CosD(bLight) * CosD(aLight)
    xLight(2) = rLight * CosD(bLight) * SinD(aLight)
    xLight(3) = rLight * SinD(bLight)

    ' Camera
    aCam = AlphaCam: bCam = BetaCam: rCam = DistCam
    xCam(1) = rCam * CosD(bCam) * CosD(aCam)
    xCam(2) = rCam * CosD(bCam) * SinD(aCam)
    xCam(3) = rCam * SinD(bCam)

    ' Screen center spatial position
    dScreen = DistScreen: aLook = 180 + aCam: bLook = -bCam
    xCenter(1) = xCam(1) + dScreen * CosD(bLook) * CosD(aLook)
    xCenter(2) = xCam(2) + dScreen * CosD(bLook) * SinD(aLook)
    xCenter(3) = xCam(3) + dScreen * SinD(bLook)

    ' Screen horizontal axis (rightwards)
    UScreen(1) = SinD(aLook): UScreen(2) = -CosD(aLook): UScreen(3) = 0
    ' Screen vertical axis (downwards)
    VScreen(1) = SinD(bLook) * CosD(aLook): VScreen(2) = SinD(bLook) * SinD(aLook): VScreen(3) = -CosD(bLook)

End Sub

Function Julia (p, q)
    ' Julia pattern
    ' p and q range from -1 to +1
    Scale = 2
    x = p * Scale: y = q * Scale
    For n% = 1 To 20
        xn = x * x - y * y + .3
        yn = 2 * x * y - .5
        If (xn * xn + yn * yn) > 4 Then Julia = 0: Exit Function
        xo = x: yo = y: x = xn: y = yn
    Next n%
    dx = xo - x: dy = yo - y
    Julia = 120 * Log((dx * dx + dy * dy) + eps) - 240
End Function

Function Light.Specular (xp(), xo())
    ' Returns the brilliance to apply to a given point on a surface
    ' The point of reflection is Xp, the normal vector to the surface is Xo
    Dim xi(3), xr(3), xv(3)
    For i = 1 To 3
        xi(i) = xp(i) - xLight(i)
    xv(i) = xCam(i) - xp(i): Next i
    Coeff = 2 * (xi(1) * xo(1) + xi(2) * xo(2) + xi(3) * xo(3))
    For i = 1 To 3: xr(i) = xi(i) - Coeff * xo(i): Next i
    NormProd = Sqr((xi(1) * xi(1) + xi(2) * xi(2) + xi(3) * xi(3)) * (xv(1) * xv(1) + xv(2) * xv(2) + xv(3) * xv(3)))
    Light.Specular = 50 + 50 * (xr(1) * xv(1) + xr(2) * xv(2) + xr(3) * xv(3)) / NormProd
End Function

Function Mandelbrot (p, q)
    ' M-Set Pattern
    ' p and q range from -1 to +1
    Scale = 4
    zx0 = p * Scale: zy0 = q * Scale
    x = zx0: y = zy0
    For NP% = 1 To 20
        xn = x * x - y * y + zx0
        yn = 2 * x * y + zy0
        If (xn * xn + yn * yn) > 4 Then Mandelbrot = 120: Exit Function
        xo = x: yo = y: x = xn: y = yn
    Next NP%
    dx = xo - x: dy = yo - y: VarMdb = dx * dx + dy * dy
    If VarMdb > 0 Then: Mandelbrot = 270 - 5 * Log(VarMdb): Else: Mandelbrot = 0
End Function

Sub Menu
    ' To allow selecting a sample of what TC-Ray can do...
    Cls
    Print "Rendering "; Name$; " took "; t!; " Seconds"
    Print
    again:
    Print "Welcome to Jark's TCRay Raytracer (ported to QB64PE by Antoni Gual)"
    Print
    Print "End of Program                    0"
    Print "Geometry                          1"
    Print "Birthday                          2"
    Print "ADN helix                         3"
    Print "Perlin Ring                       4"
    Print "Stormy Sea                        5"
    Print "Dolmen                            6"
    Print "Mandelchrist                      7"
    Print "Bike Mirrors                      8"
    Print "Water Effect                      9"
    Print "The Nephroid Cup                 10"
    Print


    Print
    Input "Select an option                : ", Choice

    Cls: 'Print "Loading space description..."
    Select Case Choice
        Case 0: Cls: EndProg
        Case 1: Init.Example1: Name$ = "Geo"
        Case 2: Init.Example2: Name$ = "Birthday"
        Case 3: Init.Example3: Name$ = "ADN"
        Case 4: Init.Example4: Name$ = "Ring"
        Case 5: Init.Example5: Name$ = "Storm"
        Case 6: Init.Example6: Name$ = "Dolmen"
        Case 7: Init.Example7: Name$ = "Christ"
        Case 8: Init.Example8: Name$ = "Bike"
        Case 9: Init.Example9: Name$ = "Water"
        Case 10: Init.Example10: Name$ = "Cup"
            '  Case 11: SelectVGA: GoTo again
        Case Else: GoTo again
    End Select


End Sub

Function Nephroid (x!, y!, Radius!)
    Nephroid = (x! ^ 2 + y! ^ 2 - 4 * (Radius! / 4) ^ 2) ^ 3 - 108 * (Radius! / 4) ^ 4 * y! ^ 2
End Function

Function Noise3! (a%, b%, C%)
    ' This controlled noise function is the basis of the 3D Perlin noise

    Noise3! = Rnd(-Rnd(-a% * b%) * Rnd(-C%) * Seed&) - .5
End Function

Function Perlin! (xp!, xq!)
    ' 2D Perlin Noise - Based on the Noise% array (shared variable)
    ' p and q range from -1 to +1

    p! = xp! / 2 + .5
    q! = xq! / 2 + .5
    If p! > 1 Then p! = 1
    If q! > 1 Then q! = 1
    If p! < 0 Then p! = 0
    If q! < 0 Then q! = 0

    ' ============================================
    ' Cosine interpolation within the noise array
    x0% = Int(p! * (Span% - 1)) + 1: y0% = Int(q! * (Span% - 1)) + 1
    dx! = p! * (Span% - 1) + 1 - x0%: dy! = q! * (Span% - 1) + 1 - y0%
    cx! = Cos(dx! * Pi / 2)
    Cy! = Cos(dy! * Pi / 2)

    If x0% = Span% And y0% = Span% Then Perlin! = Noise%(x0%, y0%) / 255: Exit Function

    If x0% = Span% And y0% < Span% Then
        z1! = Noise%(x0%, y0%)
        z2! = Noise%(x0%, y0% + 1)
        Perlin! = (z1! * Cy! + z2! * (1 - Cy!)) / 255
        Exit Function
    End If

    If y0% = Span% And x0% < Span% Then
        z1! = Noise%(x0%, y0%)
        z2! = Noise%(x0% + 1, y0%)
        Perlin! = (z1! * cx! + z2! * (1 - cx!)) / 255
        Exit Function
    End If

    z1! = Noise%(x0%, y0%)
    z2! = Noise%(x0% + 1, y0%)
    Perlin1! = z1! * cx! + z2! * (1 - cx!)
    z3! = Noise%(x0%, y0% + 1)
    z4! = Noise%(x0% + 1, y0% + 1)
    Perlin2! = z3! * cx! + z4! * (1 - cx!)
    Perlin! = (Perlin1! * Cy! + Perlin2! * (1 - Cy!)) / 255
End Function

Function Perlin3D (xt!, yt!, zt!)
    ' Perlin 3D noise. The Perlin value is calculated for every point
    ' without heavy array storage. Slow response time,
    ' but it's the only solution without EMS...

    x! = (xt! + 1) / 2 * Power2%(octaves%) + 1
    y! = (yt! + 1) / 2 * Power2%(octaves%) + 1
    z! = (zt! + 1) / 2 * Power2%(octaves%) + 1

    Rank% = 0
    ax% = 1: bx% = Span%
    ay% = 1: by% = Span%
    az% = 1: bz% = Span%
    z1% = 128: z2% = 128: z3% = 128: z4% = 128
    z5% = 128: z6% = 128: z7% = 128: z8% = 128

    Do Until (Rank% = octaves%) _OrElse (((bx% - x!) <= 1 And (x! - ax%) <= 1) And ((by% - y!) <= 1 And (y! - ay%) <= 1) And ((bz% - z!) <= 1 And (z! - az%) <= 1))
        Rank% = Rank% + 1
        Grid% = Power2%(octaves% - Rank%)

        If bx% - x! > x! - ax% Then
            bx% = bx% - Grid%
            If by% - y! > y! - ay% Then
                by% = by% - Grid%
                If bz% - z! > z! - az% Then
                    bz% = bz% - Grid%
                    z1n% = z1%
                    z2n% = (z1% + z2%) / 2
                    z3n% = (z1% + z2% + z3% + z4%) / 4
                    z4n% = (z1% + z4%) / 2
                    z5n% = (z1% + z5%) / 2
                    z6n% = (z1% + z2% + z5% + z6%) / 4
                    z7n% = (z1% + z2% + z3% + z4% + z5% + z6% + z7% + z8%) / 8
                    z8n% = (z1% + z4% + z5% + z8%) / 4
                    z1% = z1n%
                    z2% = z2n% * (1 + Noise3!(bx%, ay%, az%) * Amplitude!(Rank%))
                    z3% = z3n% * (1 + Noise3!(bx%, by%, az%) * Amplitude!(Rank%))
                    z4% = z4n% * (1 + Noise3!(ax%, by%, az%) * Amplitude!(Rank%))
                    z5% = z5n% * (1 + Noise3!(ax%, ay%, bz%) * Amplitude!(Rank%))
                    z6% = z6n% * (1 + Noise3!(bx%, ay%, bz%) * Amplitude!(Rank%))
                    z7% = z7n% * (1 + Noise3!(bx%, by%, bz%) * Amplitude!(Rank%))
                    z8% = z8n% * (1 + Noise3!(ax%, by%, bz%) * Amplitude!(Rank%))

                Else
                    az% = az% + Grid%
                    z1n% = (z1% + z5%) / 2
                    z2n% = (z1% + z2% + z5% + z6%) / 4
                    z3n% = (z1% + z2% + z3% + z4% + z5% + z6% + z7% + z8%) / 8
                    z4n% = (z1% + z4% + z5% + z8%) / 4
                    z5n% = z5%
                    z6n% = (z5% + z6%) / 2
                    z7n% = (z5% + z6% + z7% + z8%) / 4
                    z8n% = (z5% + z8%) / 2
                    z1% = z1n% * (1 + Noise3!(ax%, ay%, az%) * Amplitude!(Rank%))
                    z2% = z2n% * (1 + Noise3!(bx%, ay%, az%) * Amplitude!(Rank%))
                    z3% = z3n% * (1 + Noise3!(bx%, by%, az%) * Amplitude!(Rank%))
                    z4% = z4n% * (1 + Noise3!(ax%, by%, az%) * Amplitude!(Rank%))
                    z5% = z5n%
                    z6% = z6n% * (1 + Noise3!(bx%, ay%, bz%) * Amplitude!(Rank%))
                    z7% = z7n% * (1 + Noise3!(bx%, by%, bz%) * Amplitude!(Rank%))
                    z8% = z8n% * (1 + Noise3!(ax%, by%, bz%) * Amplitude!(Rank%))

                End If
            Else
                ay% = ay% + Grid%
                If bz% - z! > z! - az% Then
                    bz% = bz% - Grid%
                    z1n% = (z1% + z4%) / 2
                    z2n% = (z1% + z2% + z3% + z4%) / 4
                    z3n% = (z3% + z4%) / 2
                    z4n% = z4%
                    z5n% = (z1% + z4% + z5% + z8%) / 4
                    z6n% = (z1% + z2% + z3% + z4% + z5% + z6% + z7% + z8%) / 8
                    z7n% = (z4% + z3% + z7% + z8%) / 4
                    z8n% = (z4% + z8%) / 2
                    z1% = z1n% * (1 + Noise3!(ax%, ay%, az%) * Amplitude!(Rank%))
                    z2% = z2n% * (1 + Noise3!(bx%, ay%, az%) * Amplitude!(Rank%))
                    z3% = z3n% * (1 + Noise3!(bx%, by%, az%) * Amplitude!(Rank%))
                    z4% = z4n%
                    z5% = z5n% * (1 + Noise3!(ax%, ay%, bz%) * Amplitude!(Rank%))
                    z6% = z6n% * (1 + Noise3!(bx%, ay%, bz%) * Amplitude!(Rank%))
                    z7% = z7n% * (1 + Noise3!(bx%, by%, bz%) * Amplitude!(Rank%))
                    z8% = z8n% * (1 + Noise3!(ax%, by%, bz%) * Amplitude!(Rank%))
                Else
                    az% = az% + Grid%
                    z1n% = (z1% + z4% + z5% + z8%) / 4
                    z2n% = (z1% + z2% + z3% + z4% + z5% + z6% + z7% + z8%) / 8
                    z3n% = (z4% + z3% + z7% + z8%) / 4
                    z4n% = (z4% + z8%) / 2
                    z5n% = (z5% + z8%) / 2
                    z6n% = (z5% + z6% + z7% + z8%) / 4
                    z7n% = (z7% + z8%) / 2
                    z8n% = z8%
                    z1% = z1n% * (1 + Noise3!(ax%, ay%, az%) * Amplitude!(Rank%))
                    z2% = z2n% * (1 + Noise3!(bx%, ay%, az%) * Amplitude!(Rank%))
                    z3% = z3n% * (1 + Noise3!(bx%, by%, az%) * Amplitude!(Rank%))
                    z4% = z4n% * (1 + Noise3!(ax%, by%, az%) * Amplitude!(Rank%))
                    z5% = z5n% * (1 + Noise3!(ax%, ay%, bz%) * Amplitude!(Rank%))
                    z6% = z6n% * (1 + Noise3!(bx%, ay%, bz%) * Amplitude!(Rank%))
                    z7% = z7n% * (1 + Noise3!(bx%, by%, bz%) * Amplitude!(Rank%))
                    z8% = z8n%
                End If
            End If

        Else
            ax% = ax% + Grid%
            If by% - y! > y! - ay% Then
                by% = by% - Grid%
                If bz% - z! > z! - az% Then
                    bz% = bz% - Grid%
                    z1n% = (z1% + z2%) / 2
                    z2n% = z2%
                    z3n% = (z3% + z2%) / 2
                    z4n% = (z1% + z2% + z3% + z4%) / 4
                    z5n% = (z1% + z2% + z5% + z6%) / 4
                    z6n% = (z6% + z2%) / 2
                    z7n% = (z2% + z3% + z6% + z7%) / 4
                    z8n% = (z1% + z2% + z3% + z4% + z5% + z6% + z7% + z8%) / 8
                    z1% = z1n% * (1 + Noise3!(ax%, ay%, az%) * Amplitude!(Rank%))
                    z2% = z2n%
                    z3% = z3n% * (1 + Noise3!(bx%, by%, az%) * Amplitude!(Rank%))
                    z4% = z4n% * (1 + Noise3!(ax%, by%, az%) * Amplitude!(Rank%))
                    z5% = z5n% * (1 + Noise3!(ax%, ay%, bz%) * Amplitude!(Rank%))
                    z6% = z6n% * (1 + Noise3!(bx%, ay%, bz%) * Amplitude!(Rank%))
                    z7% = z7n% * (1 + Noise3!(bx%, by%, bz%) * Amplitude!(Rank%))
                    z8% = z8n% * (1 + Noise3!(ax%, by%, bz%) * Amplitude!(Rank%))
                Else
                    az% = az% + Grid%
                    z1n% = (z1% + z2% + z5% + z6%) / 4
                    z2n% = (z2% + z6%) / 2
                    z3n% = (z2% + z3% + z6% + z7%) / 4
                    z4n% = (z1% + z2% + z3% + z4% + z5% + z6% + z7% + z8%) / 8
                    z5n% = (z5% + z6%) / 2
                    z6n% = z6%
                    z7n% = (z6% + z7%) / 2
                    z8n% = (z5% + z6% + z7% + z8%) / 4
                    z1% = z1n% * (1 + Noise3!(ax%, ay%, az%) * Amplitude!(Rank%))
                    z2% = z2n% * (1 + Noise3!(bx%, ay%, az%) * Amplitude!(Rank%))
                    z3% = z3n% * (1 + Noise3!(bx%, by%, az%) * Amplitude!(Rank%))
                    z4% = z4n% * (1 + Noise3!(ax%, by%, az%) * Amplitude!(Rank%))
                    z5% = z5n% * (1 + Noise3!(ax%, ay%, bz%) * Amplitude!(Rank%))
                    z6% = z6n%
                    z7% = z7n% * (1 + Noise3!(bx%, by%, bz%) * Amplitude!(Rank%))
                    z8% = z8n% * (1 + Noise3!(ax%, by%, bz%) * Amplitude!(Rank%))
                End If
            Else
                ay% = ay% + Grid%
                If bz% - z! > z! - az% Then
                    bz% = bz% - Grid%
                    z1n% = (z1% + z2% + z3% + z4%) / 4
                    z2n% = (z3% + z2%) / 2
                    z3n% = z3%
                    z4n% = (z4% + z3%) / 2
                    z5n% = (z1% + z2% + z3% + z4% + z5% + z6% + z7% + z8%) / 8
                    z6n% = (z2% + z3% + z6% + z7%) / 4
                    z7n% = (z3% + z7%) / 2
                    z8n% = (z3% + z4% + z7% + z8%) / 4
                    z1% = z1n% * (1 + Noise3!(ax%, ay%, az%) * Amplitude!(Rank%))
                    z2% = z2n% * (1 + Noise3!(bx%, ay%, az%) * Amplitude!(Rank%))
                    z3% = z3n%
                    z4% = z4n% * (1 + Noise3!(ax%, by%, az%) * Amplitude!(Rank%))
                    z5% = z5n% * (1 + Noise3!(ax%, ay%, bz%) * Amplitude!(Rank%))
                    z6% = z6n% * (1 + Noise3!(bx%, ay%, bz%) * Amplitude!(Rank%))
                    z7% = z7n% * (1 + Noise3!(bx%, by%, bz%) * Amplitude!(Rank%))
                    z8% = z8n% * (1 + Noise3!(ax%, by%, bz%) * Amplitude!(Rank%))
                Else
                    az% = az% + Grid%
                    z1n% = (z1% + z2% + z3% + z4% + z5% + z6% + z7% + z8%) / 8
                    z2n% = (z2% + z3% + z6% + z7%) / 4
                    z3n% = (z3% + z7%) / 2
                    z4n% = (z3% + z4% + z7% + z8%) / 4
                    z5n% = (z5% + z6% + z7% + z8%) / 4
                    z6n% = (z6% + z7%) / 2
                    z7n% = z7%
                    z8n% = (z7% + z8%) / 2
                    z1% = z1n% * (1 + Noise3!(ax%, ay%, az%) * Amplitude!(Rank%))
                    z2% = z2n% * (1 + Noise3!(bx%, ay%, az%) * Amplitude!(Rank%))
                    z3% = z3n% * (1 + Noise3!(bx%, by%, az%) * Amplitude!(Rank%))
                    z4% = z4n% * (1 + Noise3!(ax%, by%, az%) * Amplitude!(Rank%))
                    z5% = z5n% * (1 + Noise3!(ax%, ay%, bz%) * Amplitude!(Rank%))
                    z6% = z6n% * (1 + Noise3!(bx%, ay%, bz%) * Amplitude!(Rank%))
                    z7% = z7n%
                    z8% = z8n% * (1 + Noise3!(ax%, by%, bz%) * Amplitude!(Rank%))
                End If
            End If

        End If
    Loop

    zp1! = z1% * (bx% - x!) + z2% * (x! - ax%)
    zp2! = z4% * (bx% - x!) + z3% * (x! - ax%)
    zp12! = zp1! * (by% - y!) + zp2! * (y! - ay%)

    zp3! = z5% * (bx% - x!) + z6% * (x! - ax%)
    zp4! = z8% * (bx% - x!) + z7% * (x! - ax%)
    zp34! = zp3! * (by% - y!) + zp4! * (y! - ay%)

    zp! = zp12! * (bz% - z!) + zp34! * (z! - az%)


    Perlin3D! = zp! / 255

End Function

Sub Pixel.Point (Nx%, Ny%, x!, y!, z!)
    Dim u(3)

    u(1) = x! - xCam(1): u(2) = y! - xCam(2): u(3) = z! - xCam(3)

    A1 = xCenter(1) - xCam(1) - xScrCenter% * UScreen(1) / ScrHeight - yScrCenter% * VScreen(1) / ScrHeight
    B1 = UScreen(1) / ScrHeight
    C1 = VScreen(1) / ScrHeight
    A2 = xCenter(2) - xCam(2) - xScrCenter% * UScreen(2) / ScrHeight - yScrCenter% * VScreen(2) / ScrHeight
    B2 = UScreen(2) / ScrHeight
    C2 = VScreen(2) / ScrHeight
    A3 = xCenter(3) - xCam(3) - xScrCenter% * UScreen(3) / ScrHeight - yScrCenter% * VScreen(3) / ScrHeight
    B3 = UScreen(3) / ScrHeight
    C3 = VScreen(3) / ScrHeight

    AlphaX1 = B1 * u(2) - B2 * u(1)
    AlphaY1 = C1 * u(2) - C2 * u(1)
    Gamma1 = -A1 * u(2) + A2 * u(1)

    AlphaX2 = B2 * u(3) - B3 * u(2)
    AlphaY2 = C2 * u(3) - C3 * u(2)
    Gamma2 = -A2 * u(3) + A3 * u(2)

    Det = AlphaX1 * AlphaY2 - AlphaX2 * AlphaY1
    Nx% = Int((Gamma1 * AlphaY2 - Gamma2 * AlphaY1) / Det + .5)
    Ny% = Int((AlphaX1 * Gamma2 - AlphaX2 * Gamma1) / Det + .5)

End Sub

Sub Probe.Ball (Id%, Ray(), Dist, xh(), xt(), xn(), FlagPlot%)
    ' Probe a plain ball with a ray
    ' Dist is the distance from the ray origin
    ' xh is the hit point in the object axis system
    ' xt is the hit point in the main axis system
    ' xn is the normal vector in the main axis system
    ' FlagPlot% is the hit indicator (0 when no intersection)

    Rx0 = Ray(1) - shape(Id%).x0
    Ry0 = Ray(2) - shape(Id%).y0
    Rz0 = Ray(3) - shape(Id%).z0
    DP = Ray(4) * Ray(4) + Ray(5) * Ray(5) + Ray(6) * Ray(6)
    DQ = 2 * (Rx0 * Ray(4) + Ry0 * Ray(5) + Rz0 * Ray(6))
    DR = Rx0 * Rx0 + Ry0 * Ry0 + Rz0 * Rz0 + Object(Id%, 35)
    Delta = DQ * DQ - 4 * DP * DR
    If Delta < 0 Then FlagPlot% = 0: Exit Sub
    Dist = -(DQ + Sqr(Delta)) * .5 / DP
    If Dist < Ray(7) Then FlagPlot% = 0: Exit Sub
    xh(1) = Rx0 + Ray(4) * Dist
    xh(2) = Ry0 + Ray(5) * Dist
    xh(3) = Rz0 + Ray(6) * Dist

    For k% = 1 To 3
        xt!(k%) = Ray(k%) + Dist * Ray(k% + 3)
        xn(k%) = xh(k%) / shape(Id%).dRef1
    Next k%
    FlagPlot% = 1
End Sub

Sub Probe.Cone (Id%, Ray(), Dist, xh(), xt(), xn(), FlagPlot%)
    ' Probe a half cone with a ray
    ' Dist is the distance from the ray origin
    ' xh is the hit point in the object axis system
    ' xt is the hit point in the main axis system
    ' xn is the normal vector in the main axis system
    ' FlagPlot% is the hit indicator (0 when no intersection)

    Dim xp!(2, 3), FlagPoint%(2), xo!(3)
    FlagPlot% = 0

    '=======================================================
    ' Ray conversion to object axis system
    Rx0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 1, 1)
    Rx0 = Rx0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 1, 2)
    Rx0 = Rx0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 1, 3)

    Ry0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 2, 1)
    Ry0 = Ry0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 2, 2)
    Ry0 = Ry0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 2, 3)

    Rz0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 3, 1)
    Rz0 = Rz0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 3, 2)
    Rz0 = Rz0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 3, 3)

    Rx = Ray(4) * Matrix(Id%, 1, 1) + Ray(5) * Matrix(Id%, 1, 2) + Ray(6) * Matrix(Id%, 1, 3)
    Ry = Ray(4) * Matrix(Id%, 2, 1) + Ray(5) * Matrix(Id%, 2, 2) + Ray(6) * Matrix(Id%, 2, 3)
    Rz = Ray(4) * Matrix(Id%, 3, 1) + Ray(5) * Matrix(Id%, 3, 2) + Ray(6) * Matrix(Id%, 3, 3)

    '=======================================================
    ' Intersection calculation in the object axis system
    DP = Rx * Rx + Ry * Ry + Object(Id%, 28) * Rz * Rz
    DQ = 2 * (Rx0 * Rx + Ry0 * Ry + Object(Id%, 28) * Rz0 * Rz)
    DR = Rx0 * Rx0 + Ry0 * Ry0 + Rz0 * Rz0 * Object(Id%, 28)
    Delta = DQ * DQ - 4 * DP * DR
    If Delta < 0 Then FlagPlot% = 0: Exit Sub
    Dist1 = -(DQ + Sqr(Delta)) * .5 / DP
    Dist2 = (Sqr(Delta) - DQ) * .5 / DP

    '=======================================================
    ' Check the intersections are not behind the screen
    FlagPoint%(1) = 1: FlagPoint%(2) = 1
    If Dist1 < Ray(7) Then FlagPoint%(1) = 0
    If Dist2 < Ray(7) Then FlagPoint%(2) = 0
    If FlagPoint%(1) = 0 And FlagPoint%(2) = 0 Then FlagPlot% = 0: Exit Sub

    xp!(1, 1) = Rx0 + Rx * Dist1: xp!(1, 2) = Ry0 + Ry * Dist1: xp!(1, 3) = Rz0 + Rz * Dist1
    xp!(2, 1) = Rx0 + Rx * Dist2: xp!(2, 2) = Ry0 + Ry * Dist2: xp!(2, 3) = Rz0 + Rz * Dist2

    '=======================================================
    ' Set limits to the shape
    For k% = 1 To 2
        If xp!(k%, 3) < 0 Or xp!(k%, 3) > shape(Id%).dRef2 Then FlagPoint%(k%) = 0
    Next k%
    If FlagPoint%(1) = 0 And FlagPoint%(2) = 0 Then FlagPlot% = 0: Exit Sub

    '=======================================================
    ' Select the intersection point (object axis system)
    If FlagPoint%(1) = 1 Then
        xh!(1) = xp!(1, 1): xh!(2) = xp!(1, 2): xh!(3) = xp!(1, 3): Dist = Dist1: FlagPlot% = 1
    Else
        xh!(1) = xp!(2, 1): xh!(2) = xp!(2, 2): xh!(3) = xp!(2, 3): Dist = Dist2: FlagPlot% = 2
    End If

    '=======================================================
    ' Calculate the normal vector
    xo!(1) = xh!(1)
    xo!(2) = xh!(2)
    xo!(3) = Object(Id%, 28) * xh!(3)
    Side = Rx * xo!(1) + Ry * xo!(2) + Rz * xo!(3)
    If Side < 0 Then Side = -1: Else: Side = 1
    ModNorm = Side * Module(xo!(1), xo!(2), xo!(3))
    If ModNorm = 0 Then Exit Sub
    For k% = 1 To 3: xo!(k%) = xo!(k%) / ModNorm: Next k%

    '=======================================================
    ' Back to main axis system
    For k% = 1 To 3
        xt!(k%) = Ray(k%) + Dist * Ray(k% + 3)
        xn!(k%) = xo!(1) * Matrix(Id%, 1, k%) + xo!(2) * Matrix(Id%, 2, k%) + xo!(3) * Matrix(Id%, 3, k%)
    Next k%


End Sub

Sub Probe.Cubic (Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%)
    ' Probe a cubic shape with a ray
    ' Dist is the distance from the ray origin
    ' xh is the hit point in the object axis system
    ' xt is the hit point in the main axis system
    ' xn is the normal vector in the main axis system
    ' FlagPlot% is the hit indicator (0 when no intersection)

    Dim xp!(3, 3), FlagPoint%(3)
    FlagPlot% = 0

    '=======================================================
    ' Ray conversion to object axis system
    Rx0# = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 1, 1)
    Rx0# = Rx0# + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 1, 2)
    Rx0# = Rx0# + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 1, 3)

    Ry0# = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 2, 1)
    Ry0# = Ry0# + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 2, 2)
    Ry0# = Ry0# + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 2, 3)

    Rz0# = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 3, 1)
    Rz0# = Rz0# + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 3, 2)
    Rz0# = Rz0# + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 3, 3)

    Rx# = Ray(4) * Matrix(Id%, 1, 1) + Ray(5) * Matrix(Id%, 1, 2) + Ray(6) * Matrix(Id%, 1, 3)
    Ry# = Ray(4) * Matrix(Id%, 2, 1) + Ray(5) * Matrix(Id%, 2, 2) + Ray(6) * Matrix(Id%, 2, 3)
    Rz# = Ray(4) * Matrix(Id%, 3, 1) + Ray(5) * Matrix(Id%, 3, 2) + Ray(6) * Matrix(Id%, 3, 3)

    '=======================================================
    ' Intersection calculation in the object axis system

    DP# = Object(Id%, 17) * Rx# * Rx# * Rx#
    DP# = DP# + Object(Id%, 18) * Ry# * Ry# * Ry#
    DP# = DP# + Object(Id%, 19) * Rz# * Rz# * Rz#
    DP# = DP# + Object(Id%, 20) * Rx# * Rx# * Ry#
    DP# = DP# + Object(Id%, 21) * Rx# * Rx# * Rz#
    DP# = DP# + Object(Id%, 23) * Ry# * Ry# * Rx#
    DP# = DP# + Object(Id%, 22) * Ry# * Ry# * Rz#
    DP# = DP# + Object(Id%, 24) * Rz# * Rz# * Rx#
    DP# = DP# + Object(Id%, 25) * Rz# * Rz# * Ry#
    DP# = DP# + Object(Id%, 16) * Rx# * Ry# * Rz#

    DQ# = DQ# + Object(Id%, 17) * 3 * Rx0# * Rx# * Rx# + Object(Id%, 18) * 3 * Ry0# * Ry# * Ry# + Object(Id%, 19) * 3 * Rz0# * Rz# * Rz#
    DQ# = DQ# + Object(Id%, 20) * Rx# * Rx# * Ry0# + Object(Id%, 21) * Rx# * Rx# * Rz0# + Object(Id%, 23) * Ry# * Ry# * Rx0#
    DQ# = DQ# + Object(Id%, 22) * Ry# * Ry# * Rz0# + Object(Id%, 24) * Rz# * Rz# * Rx0# + Object(Id%, 25) * Rz# * Rz# * Ry0#
    DQ# = DQ# + Object(Id%, 20) * 2 * Rx0# * Rx# * Ry# + Object(Id%, 21) * 2 * Rx0# * Rx# * Rz# + Object(Id%, 23) * 2 * Ry0# * Ry# * Rx#
    DQ# = DQ# + Object(Id%, 22) * 2 * Ry0# * Ry# * Rz# + Object(Id%, 24) * 2 * Rz0# * Rz# * Rx# + Object(Id%, 25) * 2 * Rz0# * Rz# * Ry#
    DQ# = DQ# + Object(Id%, 16) * Rx# * Ry# * Rz0# + Object(Id%, 16) * Rx# * Rz# * Ry0# + Object(Id%, 16) * Ry# * Rz# * Rx0#
    DQ# = DQ# + Object(Id%, 26) * Rx# * Rx# + Object(Id%, 27) * Ry# * Ry# + Object(Id%, 28) * Rz# * Rz#
    DQ# = DQ# + Object(Id%, 29) * Rx# * Ry# + Object(Id%, 30) * Ry# * Rz# + Object(Id%, 31) * Rx# * Rz#

    DR# = DR# + Object(Id%, 17) * 3 * Rx0# * Rx0# * Rx# + Object(Id%, 18) * 3 * Ry0# * Ry0# * Ry# + Object(Id%, 19) * 3 * Rz0# * Rz0# * Rz#
    DR# = DR# + Object(Id%, 20) * 2 * Rx0# * Rx# * Ry0# + Object(Id%, 21) * 2 * Rx0# * Rx# * Rz0# + Object(Id%, 23) * 2 * Rx0# * Ry0# * Ry#
    DR# = DR# + Object(Id%, 22) * 2 * Ry0# * Ry# * Rz0# + Object(Id%, 24) * 2 * Rx0# * Rz0# * Rz# + Object(Id%, 25) * 2 * Ry0# * Rz0# * Rz#
    DR# = DR# + Object(Id%, 20) * Rx0# * Rx0# * Ry# + Object(Id%, 21) * Rx0# * Rx0# * Rz# + Object(Id%, 23) * Ry0# * Ry0# * Rx#
    DR# = DR# + Object(Id%, 22) * Ry0# * Ry0# * Rz# + Object(Id%, 24) * Rz0# * Rz0# * Rx# + Object(Id%, 25) * Rz0# * Rz0# * Ry#
    DR# = DR# + Object(Id%, 16) * Rx0# * Rz0# * Ry# + Object(Id%, 16) * Ry0# * Rz0# * Rx# + Object(Id%, 16) * Rx0# * Ry0# * Rz#
    DR# = DR# + Object(Id%, 26) * 2 * Rx0# * Rx# + Object(Id%, 27) * 2 * Ry0# * Ry# + Object(Id%, 28) * 2 * Rz0# * Rz#
    DR# = DR# + Object(Id%, 29) * Ry0# * Rx# + Object(Id%, 30) * Ry0# * Rz# + Object(Id%, 31) * Rx0# * Rz# + Object(Id%, 29) * Rx0# * Ry#
    DR# = DR# + Object(Id%, 30) * Rz0# * Ry# + Object(Id%, 31) * Rz0# * Rx#
    DR# = DR# + Object(Id%, 32) * Rx# + Object(Id%, 33) * Ry# + Object(Id%, 34) * Rz#

    DS# = DS# + Object(Id%, 17) * Rx0# * Rx0# * Rx0# + Object(Id%, 18) * Ry0# * Ry0# * Ry0# + Object(Id%, 19) * Rz0# * Rz0# * Rz0#
    DS# = DS# + Object(Id%, 20) * Rx0# * Rx0# * Ry0# + Object(Id%, 21) * Rx0# * Rx0# * Rz0# + Object(Id%, 23) * Ry0# * Ry0# * Rx0#
    DS# = DS# + Object(Id%, 22) * Ry0# * Ry0# * Rz0# + Object(Id%, 24) * Rz0# * Rz0# * Rx0# + Object(Id%, 25) * Rz0# * Rz0# * Ry0#
    DS# = DS# + Object(Id%, 16) * Rx0# * Ry0# * Rz0#
    DS# = DS# + Object(Id%, 26) * Rx0# * Rx0# + Object(Id%, 27) * Ry0# * Ry0# + Object(Id%, 28) * Rz0# * Rz0#
    DS# = DS# + Object(Id%, 29) * Rx0# * Ry0# + Object(Id%, 30) * Ry0# * Rz0# + Object(Id%, 31) * Rx0# * Rz0#
    DS# = DS# + Object(Id%, 32) * Rx0# + Object(Id%, 33) * Ry0# + Object(Id%, 34) * Rz0#
    DS# = DS# + Object(Id%, 35)

    Equa3 DP#, DQ#, DR#, DS#, Dist1#, Dist2#, Dist3#, Err3%, nSol3%
    If nSol3% = 0 Then FlagPlot% = 0: Exit Sub
    If nSol3% = 1 Then Dist3# = Dist1#: Dist2# = Dist1#
    If nSol3% = 2 Then Dist3# = Dist1#
    If Dist1# > Dist2# Then Swap Dist1#, Dist2#
    If Dist1# > Dist3# Then Swap Dist1#, Dist3#
    If Dist2# > Dist3# Then Swap Dist2#, Dist3#

    '=======================================================
    ' Check the intersections are not behind the screen
    FlagPoint%(1) = 1: FlagPoint%(2) = 1: FlagPoint%(3) = 1
    If Dist1# < Ray(7) Then FlagPoint%(1) = 0
    If Dist2# < Ray(7) Then FlagPoint%(2) = 0
    If Dist3# < Ray(7) Then FlagPoint%(3) = 0
    If FlagPoint%(1) = 0 And FlagPoint%(2) = 0 And FlagPoint%(3) = 0 Then FlagPlot% = 0: Exit Sub

    xp!(1, 1) = Rx0# + Rx# * Dist1#: xp!(1, 2) = Ry0# + Ry# * Dist1#: xp!(1, 3) = Rz0# + Rz# * Dist1#
    xp!(2, 1) = Rx0# + Rx# * Dist2#: xp!(2, 2) = Ry0# + Ry# * Dist2#: xp!(2, 3) = Rz0# + Rz# * Dist2#
    xp!(3, 1) = Rx0# + Rx# * Dist3#: xp!(3, 2) = Ry0# + Ry# * Dist3#: xp!(3, 3) = Rz0# + Rz# * Dist3#

    '=======================================================
    ' Set limits to the shape
    'FOR k% = 1 TO 3
    'IF Module(xp!(k%, 1), xp!(k%, 2), xp!(k%, 3)) > Shape(Id%).dRef2 THEN FlagPoint%(k%) = 0
    'NEXT k%
    'IF FlagPoint%(1) = 0 AND FlagPoint%(2) = 0 AND FlagPoint%(3) = 0 THEN FlagPlot% = 0: EXIT SUB

    '===================================================
    ' Select the intersection point (object axis system)
    If FlagPoint%(1) = 1 Then
        xh!(1) = xp!(1, 1): xh!(2) = xp!(1, 2): xh!(3) = xp!(1, 3): Dist = Dist1#
    Else
        If FlagPoint%(2) = 1 Then
            xh!(1) = xp!(2, 1): xh!(2) = xp!(2, 2): xh!(3) = xp!(2, 3): Dist = Dist2#
        Else
            xh!(1) = xp!(3, 1): xh!(2) = xp!(3, 2): xh!(3) = xp!(3, 3): Dist = Dist3#
        End If
    End If

    '=========================
    ' Back to main axis system
    For k% = 1 To 3
        xt!(k%) = Ray(k%) + Dist * Ray(k% + 3)
    Next k%
    FlagPlot% = 1

End Sub

Sub Probe.Cylinder (Id%, Ray(), Dist, xh(), xt(), xn(), FlagPlot%)
    ' Probe a cylinder with a ray
    ' Dist is the distance from the ray origin
    ' xh is the hit point in the object axis system
    ' xt is the hit point in the main axis system
    ' xn is the normal vector in the main axis system
    ' FlagPlot% is the hit indicator (0 when no intersection)

    Dim xp!(2, 3), FlagPoint%(2), xo!(3)
    FlagPlot% = 0

    '=======================================================
    ' Ray conversion to object axis system
    Rx0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 1, 1)
    Rx0 = Rx0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 1, 2)
    Rx0 = Rx0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 1, 3)

    Ry0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 2, 1)
    Ry0 = Ry0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 2, 2)
    Ry0 = Ry0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 2, 3)

    Rz0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 3, 1)
    Rz0 = Rz0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 3, 2)
    Rz0 = Rz0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 3, 3)

    Rx = Ray(4) * Matrix(Id%, 1, 1) + Ray(5) * Matrix(Id%, 1, 2) + Ray(6) * Matrix(Id%, 1, 3)
    Ry = Ray(4) * Matrix(Id%, 2, 1) + Ray(5) * Matrix(Id%, 2, 2) + Ray(6) * Matrix(Id%, 2, 3)
    Rz = Ray(4) * Matrix(Id%, 3, 1) + Ray(5) * Matrix(Id%, 3, 2) + Ray(6) * Matrix(Id%, 3, 3)

    '=======================================================
    ' Intersection calculation in the object axis system
    DP = Rx * Rx + Ry * Ry
    DQ = 2 * (Rx0 * Rx + Ry0 * Ry)
    DR = Rx0 * Rx0 + Ry0 * Ry0 + Object(Id%, 35)

    Delta = DQ * DQ - 4 * DP * DR
    If Delta < 0 Then FlagPlot% = 0: Exit Sub
    Dist1 = -(DQ + Sqr(Delta)) * .5 / DP
    Dist2 = (Sqr(Delta) - DQ) * .5 / DP

    '=======================================================
    ' Check the intersections are not behind the screen
    FlagPoint%(1) = 1: FlagPoint%(2) = 1
    If Dist1 < Ray(7) Then FlagPoint%(1) = 0
    If Dist2 < Ray(7) Then FlagPoint%(2) = 0
    If FlagPoint%(1) = 0 And FlagPoint%(2) = 0 Then FlagPlot% = 0: Exit Sub

    xp!(1, 1) = Rx0 + Rx * Dist1: xp!(1, 2) = Ry0 + Ry * Dist1: xp!(1, 3) = Rz0 + Rz * Dist1
    xp!(2, 1) = Rx0 + Rx * Dist2: xp!(2, 2) = Ry0 + Ry * Dist2: xp!(2, 3) = Rz0 + Rz * Dist2

    '=======================================================
    ' Set limits to the shape
    For k% = 1 To 2
        If Abs(xp!(k%, 3)) > shape(Id%).dRef2 Then FlagPoint%(k%) = 0
    Next k%
    If FlagPoint%(1) = 0 And FlagPoint%(2) = 0 Then FlagPlot% = 0: Exit Sub

    '=======================================================
    ' Select the intersection point (object axis system)
    If FlagPoint%(1) = 1 Then
        xh!(1) = xp!(1, 1): xh!(2) = xp!(1, 2): xh!(3) = xp!(1, 3): Dist = Dist1: FlagPlot% = 1
    Else
        xh!(1) = xp!(2, 1): xh!(2) = xp!(2, 2): xh!(3) = xp!(2, 3): Dist = Dist2: FlagPlot% = 2
    End If

    '=======================================================
    ' Calculate the normal vector
    Side = Rx * xh!(1) + Ry * xh!(2)
    If Side < 0 Then Side = -1: Else: Side = 1
    InvModNorm = Side / shape(Id%).dRef1
    For k% = 1 To 2: xo!(k%) = xh!(k%) * InvModNorm: Next k%
    xo!(3) = 0

    '=======================================================
    ' Back to main axis system
    For k% = 1 To 3
        xt!(k%) = Ray(k%) + Dist * Ray(k% + 3)
        xn!(k%) = xo!(1) * Matrix(Id%, 1, k%) + xo!(2) * Matrix(Id%, 2, k%)
    Next k%

End Sub

Sub Probe.Floor (Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%)
    ' Probe the floor (altitude = 0) with a ray
    ' Dist is the distance from the ray origin
    ' xh is the hit point in the object axis system
    ' xt is the hit point in the main axis system
    ' xn is the normal vector in the main axis system
    ' FlagPlot% is the hit indicator (0 when no intersection)

    If Ray(6) <> 0 Then
        Dist = -Ray(3) / Ray(6)
    Else
        If Ray(3) = 0 Then: Dist = 0: Else: FlagPlot% = 0: Exit Sub
    End If
    If Dist < Ray(7) Then FlagPlot% = 0: Exit Sub
    For k% = 1 To 3
        xh(k%) = Ray(k%) + Dist * Ray(k% + 3)
        xt(k%) = xh(k%)
    Next k%
    xn(1) = 0: xn(2) = 0: xn(3) = 1
    FlagPlot% = 1
End Sub

Sub Probe.Mandelbrot (Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%)
    ' Probe a vertical M-Set with a ray
    ' Dist is the distance from the ray origin
    ' xh is the hit point in the object axis system
    ' xt is the hit point in the main axis system
    ' xn is the normal vector in the main axis system
    ' FlagPlot% is the hit indicator (0 when no intersection)

    Dim xo!(3)


    '=======================================================
    ' Ray conversion to object axis system
    Rx0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 1, 1)
    Rx0 = Rx0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 1, 2)
    Rx0 = Rx0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 1, 3)

    Ry0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 2, 1)
    Ry0 = Ry0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 2, 2)
    Ry0 = Ry0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 2, 3)

    Rz0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 3, 1)
    Rz0 = Rz0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 3, 2)
    Rz0 = Rz0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 3, 3)

    Rx = Ray(4) * Matrix(Id%, 1, 1) + Ray(5) * Matrix(Id%, 1, 2) + Ray(6) * Matrix(Id%, 1, 3)
    Ry = Ray(4) * Matrix(Id%, 2, 1) + Ray(5) * Matrix(Id%, 2, 2) + Ray(6) * Matrix(Id%, 2, 3)
    Rz = Ray(4) * Matrix(Id%, 3, 1) + Ray(5) * Matrix(Id%, 3, 2) + Ray(6) * Matrix(Id%, 3, 3)

    '=======================================================
    ' Intersection calculation in the object axis system
    DQ = Object(Id%, 32) * Rx + Object(Id%, 33) * Ry + Object(Id%, 34) * Rz
    DR = Object(Id%, 32) * Rx0 + Object(Id%, 33) * Ry0 + Object(Id%, 34) * Rz0 + Object(Id%, 35)
    If DQ <> 0 Then
        Dist = -DR / DQ
    Else
        If DR = 0 Then: Dist = 0: Else: FlagPlot% = 0: Exit Sub
    End If

    '=======================================================
    ' Check the intersection is not behind the screen
    If Dist < Ray(7) Then FlagPlot% = 0: Exit Sub

    xh!(1) = Rx0 + Rx * Dist: xh!(2) = Ry0 + Ry * Dist: xh!(3) = Rz0 + Rz * Dist: FlagPlot% = 1

    '=======================================================
    ' Set limits to the shape
    p = -xh(3): q = xh(2)
    zx0 = p: zy0 = q
    x = zx0: y = zy0
    For NP% = 1 To 20
        xn = x * x - y * y + zx0
        yn = 2 * x * y + zy0
        If (xn * xn + yn * yn) > 4 Then FlagPlot% = 0: Exit Sub
        xo = x: yo = y: x = xn: y = yn
    Next NP%
    dx = xo - x: dy = yo - y: VarMdb = dx * dx + dy * dy
    If VarMdb > 0 Then: FlagPlot% = (Int(Abs(Log(VarMdb))) Mod (4) + 1): Else: FlagPlot% = 1
    If FlagPlot% = 0 Then Exit Sub

    '=======================================================
    ' Calculate the normal vector
    xo!(1) = Object(Id%, 32): xo!(2) = Object(Id%, 33): xo!(3) = Object(Id%, 34)
    Side = Rx * xo!(1) + Ry * xo!(2) + Rz * xo!(3)
    If Side < 0 Then Side = -1: Else: Side = 1
    ModNorm = Side * Module(xo!(1), xo!(2), xo!(3))
    If ModNorm = 0 Then Exit Sub
    For k% = 1 To 3: xo!(k%) = xo!(k%) / ModNorm: Next k%

    '=======================================================
    ' Back to main axis system
    For k% = 1 To 3
        xt(k%) = Ray(k%) + Dist * Ray(k% + 3)
        xn!(k%) = xo!(1) * Matrix(Id%, 1, k%) + xo!(2) * Matrix(Id%, 2, k%) + xo!(3) * Matrix(Id%, 3, k%)
    Next k%

End Sub

Sub Probe.Object (Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%)
    ' Selects the proper probe routine
    Select Case shape(Id%).Degree
        Case 1
            Probe.Plan Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%
        Case 2
            Probe.Quadric Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%
        Case 3
            Probe.Cubic Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%
        Case 4
            Probe.Quartic Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%
        Case Else
    End Select
End Sub

Sub Probe.Plan (Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%)
    ' Probe a plan with a ray
    ' Dist is the distance from the ray origin
    ' xh is the hit point in the object axis system
    ' xt is the hit point in the main axis system
    ' xn is the normal vector in the main axis system
    ' FlagPlot% is the hit indicator (0 when no intersection)

    Dim xo!(3)

    If shape(Id%).NumType = 100 Then
        Probe.Floor Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%
        Exit Sub
    End If

    If shape(Id%).NumType = 110 Then
        Probe.Sky Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%
        Exit Sub
    End If

    If shape(Id%).NumType = 120 Then
        Probe.Mandelbrot Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%
        Exit Sub
    End If

    '=======================================================
    ' Ray conversion to object axis system
    Rx0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 1, 1)
    Rx0 = Rx0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 1, 2)
    Rx0 = Rx0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 1, 3)

    Ry0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 2, 1)
    Ry0 = Ry0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 2, 2)
    Ry0 = Ry0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 2, 3)

    Rz0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 3, 1)
    Rz0 = Rz0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 3, 2)
    Rz0 = Rz0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 3, 3)

    Rx = Ray(4) * Matrix(Id%, 1, 1) + Ray(5) * Matrix(Id%, 1, 2) + Ray(6) * Matrix(Id%, 1, 3)
    Ry = Ray(4) * Matrix(Id%, 2, 1) + Ray(5) * Matrix(Id%, 2, 2) + Ray(6) * Matrix(Id%, 2, 3)
    Rz = Ray(4) * Matrix(Id%, 3, 1) + Ray(5) * Matrix(Id%, 3, 2) + Ray(6) * Matrix(Id%, 3, 3)

    '=======================================================
    ' Intersection calculation in the object axis system
    DQ = Object(Id%, 32) * Rx + Object(Id%, 33) * Ry + Object(Id%, 34) * Rz
    DR = Object(Id%, 32) * Rx0 + Object(Id%, 33) * Ry0 + Object(Id%, 34) * Rz0 + Object(Id%, 35)
    If DQ <> 0 Then
        Dist = -DR / DQ
    Else
        If DR = 0 Then: Dist = 0: Else: FlagPlot% = 0: Exit Sub
    End If

    '=======================================================
    ' Check the intersection is not behind the screen
    If Dist < Ray(7) Then FlagPlot% = 0: Exit Sub

    xh!(1) = Rx0 + Rx * Dist: xh!(2) = Ry0 + Ry * Dist: xh!(3) = Rz0 + Rz * Dist: FlagPlot% = 1

    '=======================================================
    ' Set limits to the shape
    Select Case shape(Id%).NumType
        Case 101
            If Abs(xh!(1)) > shape(Id%).dRef1 Or Abs(xh!(2)) > shape(Id%).dRef2 Then FlagPlot% = 0
        Case 102
            If Abs(xh!(2)) > shape(Id%).dRef1 Or Abs(xh!(3)) > shape(Id%).dRef2 Then FlagPlot% = 0
        Case 103
            If Abs(xh!(1)) > shape(Id%).dRef1 Or Abs(xh!(3)) > shape(Id%).dRef2 Then FlagPlot% = 0
        Case 104
            If (xh!(1) * xh!(1) + xh!(2) * xh!(2)) > shape(Id%).dRef1 * shape(Id%).dRef1 Then FlagPlot% = 0
        Case 105
            If (xh!(3) * xh!(3) + xh!(2) * xh!(2)) > shape(Id%).dRef1 * shape(Id%).dRef1 Then FlagPlot% = 0
        Case 106
            If (xh!(1) * xh!(1) + xh!(3) * xh!(3)) > shape(Id%).dRef1 * shape(Id%).dRef1 Then FlagPlot% = 0
        Case Else
    End Select
    If FlagPlot% = 0 Then Exit Sub

    '=========================
    ' Back to main axis system
    For k% = 1 To 3
        xt(k%) = Ray(k%) + Dist * Ray(k% + 3)
    Next k%


    FlagPlot% = 1
End Sub

Sub Probe.Quadric (Id%, Ray(), Dist, xh(), xt(), xn(), FlagPlot%)
    ' Probe a quadric shape with a ray
    ' Dist is the distance from the ray origin
    ' xh is the hit point in the object axis system
    ' xt is the hit point in the main axis system
    ' xn is the normal vector in the main axis system
    ' FlagPlot% is the hit indicator (0 when no intersection)

    Dim xp!(2, 3), FlagPoint%(2), xo!(3)

    If shape(Id%).NumType = 200 Then
        Probe.Ball Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%
        Exit Sub
    End If

    If shape(Id%).NumType = 203 Then
        Probe.Cylinder Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%
        Exit Sub
    End If

    If shape(Id%).NumType = 205 Then
        Probe.Cone Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%
        Exit Sub
    End If

    'IF Shape(Id%).NumType = 210 THEN
    'Probe.Cup Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%
    'EXIT SUB
    'END IF

    FlagPlot% = 0

    '=======================================================
    ' Ray conversion to object axis system
    Rx0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 1, 1)
    Rx0 = Rx0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 1, 2)
    Rx0 = Rx0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 1, 3)

    Ry0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 2, 1)
    Ry0 = Ry0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 2, 2)
    Ry0 = Ry0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 2, 3)

    Rz0 = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 3, 1)
    Rz0 = Rz0 + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 3, 2)
    Rz0 = Rz0 + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 3, 3)

    Rx = Ray(4) * Matrix(Id%, 1, 1) + Ray(5) * Matrix(Id%, 1, 2) + Ray(6) * Matrix(Id%, 1, 3)
    Ry = Ray(4) * Matrix(Id%, 2, 1) + Ray(5) * Matrix(Id%, 2, 2) + Ray(6) * Matrix(Id%, 2, 3)
    Rz = Ray(4) * Matrix(Id%, 3, 1) + Ray(5) * Matrix(Id%, 3, 2) + Ray(6) * Matrix(Id%, 3, 3)

    '=======================================================
    ' Intersection calculation in the object axis system
    DP = Object(Id%, 26) * Rx * Rx + Object(Id%, 27) * Ry * Ry + Object(Id%, 28) * Rz * Rz
    DP = DP + Object(Id%, 29) * Rx * Ry + Object(Id%, 30) * Ry * Rz + Object(Id%, 31) * Rx * Rz
    DQ = 2 * (Object(Id%, 26) * Rx0 * Rx + Object(Id%, 27) * Ry0 * Ry + Object(Id%, 28) * Rz0 * Rz)
    DQ = DQ + Object(Id%, 29) * (Rx0 * Ry + Rx * Ry0) + Object(Id%, 30) * (Ry0 * Rz + Ry * Rz0) + Object(Id%, 31) * (Rx0 * Rz + Rx * Rz0)
    DQ = DQ + Object(Id%, 32) * Rx + Object(Id%, 33) * Ry + Object(Id%, 34) * Rz
    DR = Rx0 * Rx0 * Object(Id%, 26) + Ry0 * Ry0 * Object(Id%, 27) + Rz0 * Rz0 * Object(Id%, 28) + Object(Id%, 29) * Rx0 * Ry0 + Object(Id%, 30) * Ry0 * Rz0 + Object(Id%, 31) * Rx0 * Rz0
    DR = DR + Object(Id%, 32) * Rx0 + Object(Id%, 33) * Ry0 + Object(Id%, 34) * Rz0 + Object(Id%, 35)

    If DP <> 0 Then
        Delta = DQ * DQ - 4 * DP * DR
        If Delta >= 0 Then
            Dist1 = -(DQ + Sqr(Delta)) * .5 / DP
            Dist2 = (Sqr(Delta) - DQ) * .5 / DP
        Else
            FlagPlot% = 0: Exit Sub
        End If
    Else
        If DQ <> 0 Then
            Dist1 = -DR / DQ: Dist2 = Dist1
        Else
            If DR = 0 Then: Dist1 = 0: Dist2 = 0: Else: FlagPlot% = 0: Exit Sub
        End If
    End If

    '=======================================================
    ' Check the intersections are not behind the screen
    FlagPoint%(1) = 1: FlagPoint%(2) = 1
    If Dist1 < Ray(7) Then FlagPoint%(1) = 0
    If Dist2 < Ray(7) Then FlagPoint%(2) = 0
    If FlagPoint%(1) = 0 And FlagPoint%(2) = 0 Then FlagPlot% = 0: Exit Sub

    xp!(1, 1) = Rx0 + Rx * Dist1: xp!(1, 2) = Ry0 + Ry * Dist1: xp!(1, 3) = Rz0 + Rz * Dist1
    xp!(2, 1) = Rx0 + Rx * Dist2: xp!(2, 2) = Ry0 + Ry * Dist2: xp!(2, 3) = Rz0 + Rz * Dist2

    '=======================================================
    ' Set limits to the shape
    For k% = 1 To 2
        Select Case shape(Id%).NumType
            Case 201
                If Abs(xp!(k%, 1)) > shape(Id%).dRef2 Or Abs(xp!(k%, 2)) > shape(Id%).dRef2 Or Abs(xp!(k%, 3)) > shape(Id%).dRef2 Then FlagPoint%(k%) = 0
            Case 202
            Case 206
                If (xp!(k%, 1) ^ 2 + xp!(k%, 2) ^ 2) > shape(Id%).dRef1 ^ 2 Then FlagPoint%(k%) = 0
            Case 210
                If xp!(k%, 3) > 0 Then FlagPoint%(k%) = 0
            Case Else
                If Abs(xp!(k%, 3)) > shape(Id%).dRef2 Then FlagPoint%(k%) = 0
        End Select
    Next k%
    If FlagPoint%(1) = 0 And FlagPoint%(2) = 0 Then FlagPlot% = 0: Exit Sub

    '===================================================
    ' Select the intersection point (object axis system)
    If FlagPoint%(1) = 1 Then
        xh!(1) = xp!(1, 1): xh!(2) = xp!(1, 2): xh!(3) = xp!(1, 3): Dist = Dist1: FlagPlot% = 1
    Else
        xh!(1) = xp!(2, 1): xh!(2) = xp!(2, 2): xh!(3) = xp!(2, 3): Dist = Dist2: FlagPlot% = 2
    End If

    '=========================
    ' Back to main axis system
    For k% = 1 To 3
        xt(k%) = Ray(k%) + Dist * Ray(k% + 3)
    Next k%
End Sub

Sub Probe.Quartic (Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%)
    ' Probe a quartic shape with a ray
    ' Dist is the distance from the ray origin
    ' xh is the hit point in the object axis system
    ' xt is the hit point in the main axis system
    ' xn is the normal vector in the main axis system
    ' FlagPlot% is the hit indicator (0 when no intersection)

    Dim xp!(4, 3), FlagPoint%(4), xo!(3)
    FlagPlot% = 0

    '=======================================================
    ' Ray conversion to object axis system
    Rx0# = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 1, 1)
    Rx0# = Rx0# + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 1, 2)
    Rx0# = Rx0# + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 1, 3)

    Ry0# = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 2, 1)
    Ry0# = Ry0# + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 2, 2)
    Ry0# = Ry0# + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 2, 3)

    Rz0# = (Ray(1) - shape(Id%).x0) * Matrix(Id%, 3, 1)
    Rz0# = Rz0# + (Ray(2) - shape(Id%).y0) * Matrix(Id%, 3, 2)
    Rz0# = Rz0# + (Ray(3) - shape(Id%).z0) * Matrix(Id%, 3, 3)

    Rx# = Ray(4) * Matrix(Id%, 1, 1) + Ray(5) * Matrix(Id%, 1, 2) + Ray(6) * Matrix(Id%, 1, 3)
    Ry# = Ray(4) * Matrix(Id%, 2, 1) + Ray(5) * Matrix(Id%, 2, 2) + Ray(6) * Matrix(Id%, 2, 3)
    Rz# = Ray(4) * Matrix(Id%, 3, 1) + Ray(5) * Matrix(Id%, 3, 2) + Ray(6) * Matrix(Id%, 3, 3)

    '=======================================================
    ' Intersection calculation in the object axis system


    DP# = Object(Id%, 1) * Rx# * Rx# * Rx# * Rx#
    DP# = DP# + Object(Id%, 2) * Ry# * Ry# * Ry# * Ry#
    DP# = DP# + Object(Id%, 3) * Rz# * Rz# * Rz# * Rz#
    DP# = DP# + Object(Id%, 13) * Rx# * Rx# * Ry# * Rz#
    DP# = DP# + Object(Id%, 14) * Ry# * Ry# * Rx# * Rz#
    DP# = DP# + Object(Id%, 15) * Rz# * Rz# * Rx# * Ry#
    DP# = DP# + Object(Id%, 10) * Rx# * Rx# * Ry# * Ry#
    DP# = DP# + Object(Id%, 11) * Rx# * Rx# * Rz# * Rz#
    DP# = DP# + Object(Id%, 12) * Ry# * Ry# * Rz# * Rz#
    DP# = DP# + Object(Id%, 4) * Rx# * Rx# * Rx# * Ry# + Object(Id%, 5) * Rx# * Rx# * Rx# * Rz#
    DP# = DP# + Object(Id%, 6) * Ry# * Ry# * Ry# * Rz# + Object(Id%, 7) * Ry# * Ry# * Ry# * Rx#
    DP# = DP# + Object(Id%, 8) * Rz# * Rz# * Rz# * Rx# + Object(Id%, 9) * Rz# * Rz# * Rz# * Ry#

    DQ# = Object(Id%, 1) * 4 * Rx0# * Rx# * Rx# * Rx#
    DQ# = DQ# + Object(Id%, 2) * 4 * Ry0# * Ry# * Ry# * Ry#
    DQ# = DQ# + Object(Id%, 3) * 4 * Rz0# * Rz# * Rz# * Rz#
    DQ# = DQ# + Object(Id%, 13) * Rx# * Rx# * Ry# * Rz0#
    DQ# = DQ# + Object(Id%, 14) * Ry# * Ry# * Rx# * Rz0#
    DQ# = DQ# + Object(Id%, 15) * Rz# * Rz# * Rx# * Ry0#
    DQ# = DQ# + Object(Id%, 13) * Rx# * Rx# * Ry0# * Rz#
    DQ# = DQ# + Object(Id%, 14) * Ry# * Ry# * Rx0# * Rz#
    DQ# = DQ# + Object(Id%, 15) * Rz# * Rz# * Rx0# * Ry#
    DQ# = DQ# + Object(Id%, 13) * 2 * Rx0# * Rx# * Ry# * Rz#
    DQ# = DQ# + Object(Id%, 14) * 2 * Ry0# * Ry# * Rx# * Rz#
    DQ# = DQ# + Object(Id%, 15) * 2 * Rz0# * Rz# * Rx# * Ry#
    DQ# = DQ# + Object(Id%, 10) * 2 * Rx# * Rx# * Ry0# * Ry#
    DQ# = DQ# + Object(Id%, 11) * 2 * Rx# * Rx# * Rz0# * Rz#
    DQ# = DQ# + Object(Id%, 12) * 2 * Ry# * Ry# * Rz0# * Rz#
    DQ# = DQ# + Object(Id%, 10) * 2 * Ry# * Ry# * Rx0# * Rx#
    DQ# = DQ# + Object(Id%, 11) * 2 * Rz# * Rz# * Rx0# * Rx#
    DQ# = DQ# + Object(Id%, 12) * 2 * Rz# * Rz# * Ry0# * Ry#
    DQ# = DQ# + Object(Id%, 4) * Rx# * Rx# * Rx# * Ry0#
    DQ# = DQ# + Object(Id%, 5) * Rx# * Rx# * Rx# * Rz0#
    DQ# = DQ# + Object(Id%, 6) * Ry# * Ry# * Ry# * Rz0#
    DQ# = DQ# + Object(Id%, 7) * Ry# * Ry# * Ry# * Rx0#
    DQ# = DQ# + Object(Id%, 9) * Rz# * Rz# * Rz# * Rx0#
    DQ# = DQ# + Object(Id%, 9) * Rz# * Rz# * Rz# * Ry0#
    DQ# = DQ# + Object(Id%, 4) * 3 * Rx0# * Rx# * Rx# * Ry#
    DQ# = DQ# + Object(Id%, 5) * 3 * Rx0# * Rx# * Rx# * Rz#
    DQ# = DQ# + Object(Id%, 6) * 3 * Ry0# * Ry# * Ry# * Rz#
    DQ# = DQ# + Object(Id%, 7) * 3 * Ry0# * Ry# * Ry# * Rx#
    DQ# = DQ# + Object(Id%, 8) * 3 * Rz0# * Rz# * Rz# * Rx#
    DQ# = DQ# + Object(Id%, 9) * 3 * Rz0# * Rz# * Rz# * Ry#
    DQ# = DQ# + Object(Id%, 17) * Rx# * Rx# * Rx#
    DQ# = DQ# + Object(Id%, 18) * Ry# * Ry# * Ry#
    DQ# = DQ# + Object(Id%, 19) * Rz# * Rz# * Rz#
    DQ# = DQ# + Object(Id%, 20) * Rx# * Rx# * Ry#
    DQ# = DQ# + Object(Id%, 21) * Rx# * Rx# * Rz#
    DQ# = DQ# + Object(Id%, 23) * Ry# * Ry# * Rx#
    DQ# = DQ# + Object(Id%, 22) * Ry# * Ry# * Rz#
    DQ# = DQ# + Object(Id%, 24) * Rz# * Rz# * Rx#
    DQ# = DQ# + Object(Id%, 25) * Rz# * Rz# * Ry#
    DQ# = DQ# + Object(Id%, 16) * Rx# * Ry# * Rz#

    DR# = Object(Id%, 1) * 6 * Rx0# * Rx0# * Rx# * Rx#
    DR# = DR# + Object(Id%, 2) * 6 * Ry0# * Ry0# * Ry# * Ry#
    DR# = DR# + Object(Id%, 3) * 6 * Rz0# * Rz0# * Rz# * Rz#
    DR# = DR# + Object(Id%, 13) * 2 * Rx0# * Rx# * Ry# * Rz0#
    DR# = DR# + Object(Id%, 14) * 2 * Ry0# * Ry# * Rx# * Rz0#
    DR# = DR# + Object(Id%, 15) * 2 * Rz0# * Rz# * Rx# * Ry0#
    DR# = DR# + Object(Id%, 13) * 2 * Rx0# * Rx# * Ry0# * Rz#
    DR# = DR# + Object(Id%, 14) * 2 * Ry0# * Ry# * Rx0# * Rz#
    DR# = DR# + Object(Id%, 15) * 2 * Rz0# * Rz# * Rx0# * Ry#
    DR# = DR# + Object(Id%, 13) * Rx# * Rx# * Ry0# * Rz0#
    DR# = DR# + Object(Id%, 14) * Ry# * Ry# * Rx0# * Rz0#
    DR# = DR# + Object(Id%, 15) * Rz# * Rz# * Rx0# * Ry0#
    DR# = DR# + Object(Id%, 13) * Rx0# * Rx0# * Ry# * Rz#
    DR# = DR# + Object(Id%, 14) * Ry0# * Ry0# * Rx# * Rz#
    DR# = DR# + Object(Id%, 15) * Rz0# * Rz0# * Rx# * Ry#
    DR# = DR# + Object(Id%, 10) * Rx# * Rx# * Ry0# * Ry0#
    DR# = DR# + Object(Id%, 11) * Rx# * Rx# * Rz0# * Rz0#
    DR# = DR# + Object(Id%, 12) * Ry# * Ry# * Rz0# * Rz0#
    DR# = DR# + Object(Id%, 10) * Rx0# * Rx0# * Ry# * Ry#
    DR# = DR# + Object(Id%, 11) * Rx0# * Rx0# * Rz# * Rz#
    DR# = DR# + Object(Id%, 12) * Ry0# * Ry0# * Rz# * Rz#
    DR# = DR# + Object(Id%, 10) * 4 * Rx0# * Rx# * Ry0# * Ry#
    DR# = DR# + Object(Id%, 11) * 4 * Rx0# * Rx# * Rz0# * Rz#
    DR# = DR# + Object(Id%, 12) * 4 * Ry0# * Ry# * Rz0# * Rz#
    DR# = DR# + Object(Id%, 4) * 3 * Rx0# * Rx# * Rx# * Ry0#
    DR# = DR# + Object(Id%, 5) * 3 * Rx0# * Rx# * Rx# * Rz0#
    DR# = DR# + Object(Id%, 6) * 3 * Ry0# * Ry# * Ry# * Rz0#
    DR# = DR# + Object(Id%, 7) * 3 * Ry0# * Ry# * Ry# * Rx0#
    DR# = DR# + Object(Id%, 8) * 3 * Rz0# * Rz# * Rz# * Rx0#
    DR# = DR# + Object(Id%, 9) * 3 * Rz0# * Rz# * Rz# * Ry0#
    DR# = DR# + Object(Id%, 4) * 3 * Rx0# * Rx0# * Rx# * Ry#
    DR# = DR# + Object(Id%, 5) * 3 * Rx0# * Rx0# * Rx# * Rz#
    DR# = DR# + Object(Id%, 6) * 3 * Ry0# * Ry0# * Ry# * Rz#
    DR# = DR# + Object(Id%, 7) * 3 * Ry0# * Ry0# * Ry# * Rx#
    DR# = DR# + Object(Id%, 8) * 3 * Rz0# * Rz0# * Rz# * Rx#
    DR# = DR# + Object(Id%, 9) * 3 * Rz0# * Rz0# * Rz# * Ry#
    DR# = DR# + Object(Id%, 17) * 3 * Rx0# * Rx# * Rx#
    DR# = DR# + Object(Id%, 18) * 3 * Ry0# * Ry# * Ry#
    DR# = DR# + Object(Id%, 19) * 3 * Rz0# * Rz# * Rz#
    DR# = DR# + Object(Id%, 20) * Rx# * Rx# * Ry0#
    DR# = DR# + Object(Id%, 21) * Rx# * Rx# * Rz0#
    DR# = DR# + Object(Id%, 23) * Ry# * Ry# * Rx0#
    DR# = DR# + Object(Id%, 22) * Ry# * Ry# * Rz0#
    DR# = DR# + Object(Id%, 24) * Rz# * Rz# * Rx0#
    DR# = DR# + Object(Id%, 25) * Rz# * Rz# * Ry0#
    DR# = DR# + Object(Id%, 20) * 2 * Rx0# * Rx# * Ry#
    DR# = DR# + Object(Id%, 21) * 2 * Rx0# * Rx# * Rz#
    DR# = DR# + Object(Id%, 23) * 2 * Ry0# * Ry# * Rx#
    DR# = DR# + Object(Id%, 22) * 2 * Ry0# * Ry# * Rz#
    DR# = DR# + Object(Id%, 24) * 2 * Rz0# * Rz# * Rx#
    DR# = DR# + Object(Id%, 25) * 2 * Rz0# * Rz# * Ry#
    DR# = DR# + Object(Id%, 16) * Rx# * Ry# * Rz0#
    DR# = DR# + Object(Id%, 16) * Rx# * Rz# * Ry0#
    DR# = DR# + Object(Id%, 16) * Ry# * Rz# * Rx0#
    DR# = DR# + Object(Id%, 26) * Rx# * Rx#
    DR# = DR# + Object(Id%, 27) * Ry# * Ry#
    DR# = DR# + Object(Id%, 28) * Rz# * Rz#
    DR# = DR# + Object(Id%, 29) * Rx# * Ry#
    DR# = DR# + Object(Id%, 30) * Ry# * Rz#
    DR# = DR# + Object(Id%, 31) * Rx# * Rz#

    DS# = Object(Id%, 1) * 4 * Rx0# * Rx0# * Rx0# * Rx#
    DS# = DS# + Object(Id%, 2) * 4 * Ry0# * Ry0# * Ry0# * Ry#
    DS# = DS# + Object(Id%, 3) * 4 * Rz0# * Rz0# * Rz0# * Rz#
    DS# = DS# + Object(Id%, 13) * 2 * Rx0# * Rx# * Ry0# * Rz0#
    DS# = DS# + Object(Id%, 14) * 2 * Ry0# * Ry# * Rx0# * Rz0#
    DS# = DS# + Object(Id%, 15) * 2 * Rz0# * Rz# * Rx0# * Ry0#
    DS# = DS# + Object(Id%, 13) * Rx0# * Rx0# * Ry# * Rz0#
    DS# = DS# + Object(Id%, 14) * 2 * Ry0# * Ry0# * Rx# * Rz0#
    DS# = DS# + Object(Id%, 15) * 2 * Rz0# * Rz0# * Rx# * Ry0#
    DS# = DS# + Object(Id%, 13) * Rx0# * Rx0# * Ry0# * Rz#
    DS# = DS# + Object(Id%, 14) * 2 * Ry0# * Ry0# * Rx0# * Rz#
    DS# = DS# + Object(Id%, 15) * 2 * Rz0# * Rz0# * Rx0# * Ry#
    DS# = DS# + Object(Id%, 10) * 2 * Rx0# * Rx# * Ry0# * Ry0#
    DS# = DS# + Object(Id%, 11) * 2 * Rx0# * Rx# * Rz0# * Rz0#
    DS# = DS# + Object(Id%, 12) * 2 * Ry0# * Ry# * Rz0# * Rz0#
    DS# = DS# + Object(Id%, 10) * 2 * Ry0# * Ry# * Rx0# * Rx0#
    DS# = DS# + Object(Id%, 11) * 2 * Rz0# * Rz# * Rx0# * Rx0#
    DS# = DS# + Object(Id%, 12) * 2 * Rz0# * Rz# * Ry0# * Ry0#
    DS# = DS# + Object(Id%, 4) * 3 * Rx0# * Rx0# * Rx# * Ry0#
    DS# = DS# + Object(Id%, 5) * 3 * Rx0# * Rx0# * Rx# * Rz0#
    DS# = DS# + Object(Id%, 6) * 3 * Ry0# * Ry0# * Ry# * Rz0#
    DS# = DS# + Object(Id%, 7) * 3 * Ry0# * Ry0# * Rx0# * Ry#
    DS# = DS# + Object(Id%, 8) * 3 * Rz0# * Rz0# * Rx0# * Rz#
    DS# = DS# + Object(Id%, 9) * 3 * Rz0# * Rz0# * Ry0# * Rz#
    DS# = DS# + Object(Id%, 4) * Rx0# * Rx0# * Rx0# * Ry#
    DS# = DS# + Object(Id%, 5) * Rx0# * Rx0# * Rx0# * Rz#
    DS# = DS# + Object(Id%, 6) * Ry0# * Ry0# * Ry0# * Rz#
    DS# = DS# + Object(Id%, 7) * Ry0# * Ry0# * Ry0# * Rx#
    DS# = DS# + Object(Id%, 8) * Rz0# * Rz0# * Rz0# * Rx#
    DS# = DS# + Object(Id%, 9) * Rz0# * Rz0# * Rz0# * Ry#
    DS# = DS# + Object(Id%, 17) * 3 * Rx0# * Rx0# * Rx#
    DS# = DS# + Object(Id%, 18) * 3 * Ry0# * Ry0# * Ry#
    DS# = DS# + Object(Id%, 19) * 3 * Rz0# * Rz0# * Rz#
    DS# = DS# + Object(Id%, 20) * 2 * Rx0# * Rx# * Ry0#
    DS# = DS# + Object(Id%, 21) * 2 * Rx0# * Rx# * Rz0#
    DS# = DS# + Object(Id%, 23) * 2 * Rx0# * Ry0# * Ry#
    DS# = DS# + Object(Id%, 22) * 2 * Ry0# * Ry# * Rz0#
    DS# = DS# + Object(Id%, 24) * 2 * Rx0# * Rz0# * Rz#
    DS# = DS# + Object(Id%, 25) * 2 * Ry0# * Rz0# * Rz#
    DS# = DS# + Object(Id%, 20) * Rx0# * Rx0# * Ry#
    DS# = DS# + Object(Id%, 21) * Rx0# * Rx0# * Rz#
    DS# = DS# + Object(Id%, 23) * Ry0# * Ry0# * Rx#
    DS# = DS# + Object(Id%, 22) * Ry0# * Ry0# * Rz#
    DS# = DS# + Object(Id%, 24) * Rz0# * Rz0# * Rx#
    DS# = DS# + Object(Id%, 25) * Rz0# * Rz0# * Ry#
    DS# = DS# + Object(Id%, 16) * Rx0# * Rz0# * Ry#
    DS# = DS# + Object(Id%, 16) * Ry0# * Rz0# * Rx#
    DS# = DS# + Object(Id%, 16) * Rx0# * Ry0# * Rz#
    DS# = DS# + Object(Id%, 26) * 2 * Rx0# * Rx#
    DS# = DS# + Object(Id%, 27) * 2 * Ry0# * Ry#
    DS# = DS# + Object(Id%, 28) * 2 * Rz0# * Rz#
    DS# = DS# + Object(Id%, 29) * Ry0# * Rx#
    DS# = DS# + Object(Id%, 30) * Ry0# * Rz#
    DS# = DS# + Object(Id%, 31) * Rx0# * Rz#
    DS# = DS# + Object(Id%, 29) * Rx0# * Ry#
    DS# = DS# + Object(Id%, 30) * Rz0# * Ry#
    DS# = DS# + Object(Id%, 31) * Rz0# * Rx#
    DS# = DS# + Object(Id%, 32) * Rx#
    DS# = DS# + Object(Id%, 33) * Ry#
    DS# = DS# + Object(Id%, 34) * Rz#

    dt# = Object(Id%, 1) * Rx0# * Rx0# * Rx0# * Rx0#
    dt# = dt# + Object(Id%, 2) * Ry0# * Ry0# * Ry0# * Ry0#
    dt# = dt# + Object(Id%, 3) * Rz0# * Rz0# * Rz0# * Rz0#
    dt# = dt# + Object(Id%, 13) * Rx0# * Rx0# * Ry0# * Rz0#
    dt# = dt# + Object(Id%, 14) * Ry0# * Ry0# * Rx0# * Rz0#
    dt# = dt# + Object(Id%, 15) * Rz0# * Rz0# * Rx0# * Ry0#
    dt# = dt# + Object(Id%, 10) * Rx0# * Rx0# * Ry0# * Ry0#
    dt# = dt# + Object(Id%, 11) * Rx0# * Rx0# * Rz0# * Rz0#
    dt# = dt# + Object(Id%, 12) * Ry0# * Ry0# * Rz0# * Rz0#
    dt# = dt# + Object(Id%, 4) * Rx0# * Rx0# * Rx0# * Ry0#
    dt# = dt# + Object(Id%, 5) * Rx0# * Rx0# * Rx0# * Rz0#
    dt# = dt# + Object(Id%, 6) * Ry0# * Ry0# * Ry0# * Rz0#
    dt# = dt# + Object(Id%, 7) * Ry0# * Ry0# * Ry0# * Rx0#
    dt# = dt# + Object(Id%, 8) * Rz0# * Rz0# * Rz0# * Rx0#
    dt# = dt# + Object(Id%, 9) * Rz0# * Rz0# * Rz0# * Ry0#
    dt# = dt# + Object(Id%, 17) * Rx0# * Rx0# * Rx0#
    dt# = dt# + Object(Id%, 18) * Ry0# * Ry0# * Ry0#
    dt# = dt# + Object(Id%, 19) * Rz0# * Rz0# * Rz0#
    dt# = dt# + Object(Id%, 20) * Rx0# * Rx0# * Ry0#
    dt# = dt# + Object(Id%, 21) * Rx0# * Rx0# * Rz0#
    dt# = dt# + Object(Id%, 23) * Ry0# * Ry0# * Rx0#
    dt# = dt# + Object(Id%, 22) * Ry0# * Ry0# * Rz0#
    dt# = dt# + Object(Id%, 24) * Rz0# * Rz0# * Rx0#
    dt# = dt# + Object(Id%, 25) * Rz0# * Rz0# * Ry0#
    dt# = dt# + Object(Id%, 16) * Rx0# * Ry0# * Rz0#
    dt# = dt# + Object(Id%, 26) * Rx0# * Rx0#
    dt# = dt# + Object(Id%, 27) * Ry0# * Ry0#
    dt# = dt# + Object(Id%, 28) * Rz0# * Rz0#
    dt# = dt# + Object(Id%, 29) * Rx0# * Ry0#
    dt# = dt# + Object(Id%, 30) * Ry0# * Rz0#
    dt# = dt# + Object(Id%, 31) * Rx0# * Rz0#
    dt# = dt# + Object(Id%, 32) * Rx0#
    dt# = dt# + Object(Id%, 33) * Ry0#
    dt# = dt# + Object(Id%, 34) * Rz0#
    dt# = dt# + Object(Id%, 35)

    Equa4 DP#, DQ#, DR#, DS#, dt#, Dist1#, Dist2#, Dist3#, Dist4#, Err4%, nSol4%

    If nSol4% = 0 Then FlagPlot% = 0: Exit Sub
    If nSol4% = 1 Then Dist4# = Dist1#: Dist3# = Dist1#: Dist2# = Dist1#
    If nSol4% = 2 Then Dist4# = Dist2#: Dist3# = Dist1#
    If nSol4% = 3 Then Dist4# = Dist3#
    If Dist1# > Dist2# Then Swap Dist1#, Dist2#
    If Dist1# > Dist3# Then Swap Dist1#, Dist3#
    If Dist1# > Dist4# Then Swap Dist1#, Dist4#
    If Dist2# > Dist3# Then Swap Dist2#, Dist3#
    If Dist2# > Dist4# Then Swap Dist2#, Dist4#
    If Dist3# > Dist4# Then Swap Dist3#, Dist4#

    '=======================================================
    ' Check the intersections are not behind the screen
    FlagPoint%(1) = 1: FlagPoint%(2) = 1: FlagPoint%(3) = 1: FlagPoint%(4) = 1
    If Dist1# < Ray(7) Then FlagPoint%(1) = 0
    If Dist2# < Ray(7) Then FlagPoint%(2) = 0
    If Dist3# < Ray(7) Then FlagPoint%(3) = 0
    If Dist4# < Ray(7) Then FlagPoint%(4) = 0

    If (FlagPoint%(1) = 0 And FlagPoint%(2) = 0 And FlagPoint%(3) = 0 And FlagPoint%(4) = 0) Then FlagPlot% = 0: Exit Sub
    If FlagPoint%(1) = 0 Then FlagPlot% = 0: Exit Sub

    xp!(1, 1) = Rx0# + Rx# * Dist1#: xp!(1, 2) = Ry0# + Ry# * Dist1#: xp!(1, 3) = Rz0# + Rz# * Dist1#
    xp!(2, 1) = Rx0# + Rx# * Dist2#: xp!(2, 2) = Ry0# + Ry# * Dist2#: xp!(2, 3) = Rz0# + Rz# * Dist2#
    xp!(3, 1) = Rx0# + Rx# * Dist3#: xp!(3, 2) = Ry0# + Ry# * Dist3#: xp!(3, 3) = Rz0# + Rz# * Dist3#
    xp!(4, 1) = Rx0# + Rx# * Dist4#: xp!(4, 2) = Ry0# + Ry# * Dist4#: xp!(4, 3) = Rz0# + Rz# * Dist4#

    '=======================================================
    ' Select the intersection point (object axis system)
    If FlagPoint%(1) = 1 Then
        xh!(1) = xp!(1, 1): xh!(2) = xp!(1, 2): xh!(3) = xp!(1, 3): Dist = Dist1#
    Else
        If FlagPoint%(2) = 1 Then
            xh!(1) = xp!(2, 1): xh!(2) = xp!(2, 2): xh!(3) = xp!(2, 3): Dist = Dist2#
        Else
            If FlagPoint%(3) = 1 Then
                xh!(1) = xp!(3, 1): xh!(2) = xp!(3, 2): xh!(3) = xp!(3, 3): Dist = Dist3#
            Else
                xh!(1) = xp!(4, 1): xh!(2) = xp!(4, 2): xh!(3) = xp!(4, 3): Dist = Dist4#
            End If
        End If
    End If

    '=========================
    ' Back to main axis system
    For k% = 1 To 3
        xt(k%) = Ray(k%) + Dist * Ray(k% + 3)
    Next k%
    FlagPlot% = 1
End Sub

Sub Probe.Ray (Ray(), ShapeHit%, DistHit, xhHit(), xtHit(), xnHit(), FlagHit%)
    ' Probe space with a ray
    ' DistHit is the distance to the closest object
    ' xhHit is the hit point in the object axis system
    ' xtHit is the hit point in the main axis system
    ' xnHit is the normal vector in the main axis system
    ' FlagHit% is the hit indicator (0 when no intersection)

    Dim xh(3), xt(3), xn(3), RayHit(7)
    FlagHit% = 0 ' Reset hit flag
    DistHit = 1E+12 ' Start from faraway

    ' Probe all objects with the ray
    For Id% = 1 To NbObjects%
        Probe.Object Id%, Ray(), Dist, xh(), xt(), xn(), FlagPlot%
        If FlagPlot% <> 0 And Dist < DistHit And Dist >= 0 Then
            FlagHit% = FlagPlot%: ShapeHit% = Id%: DistHit = Dist
            For k% = 1 To 3: xhHit(k%) = xh(k%): xtHit(k%) = xt(k%): xnHit(k%) = xn(k%): Next k%
            For k% = 1 To 7: RayHit(k%) = Ray(k%): Next k%
        End If
    Next Id%

    ' Quit if nothing was hit
    If FlagHit% = 0 Then Exit Sub

    ' Get the normal vector
    Select Case shape(ShapeHit%).NumType
        Case 100, 110, 120, 200, 203, 205
            ' Optimised routines : don't recalculate normal
        Case Else
            ' Get the normal vector for the rest
            Calc.Normal ShapeHit%, xhHit(), xnHit(), RayHit()
    End Select

End Sub

Sub Probe.Sky (Id%, Ray(), Dist, xh!(), xt!(), xn!(), FlagPlot%)
    ' Probe the sky (Altitude = z0) with a ray
    ' Dist is the distance from the ray origin
    ' xh is the hit point in the object axis system
    ' xt is the hit point in the main axis system
    ' xn is the normal vector in the main axis system
    ' FlagPlot% is the hit indicator (0 when no intersection)

    'Rx0 = Ray(1): Ry0 = Ray(2):
    Rz0 = Ray(3) - shape(Id%).z0

    If Ray(6) <> 0 Then
        Dist = -Rz0 / Ray(6)
    Else
        If Ray(3) = 0 Then: Dist = 0: Else: FlagPlot% = 0: Exit Sub
    End If
    If Dist < Ray(7) Then FlagPlot% = 0: Exit Sub
    For k% = 1 To 3
        xh(k%) = Ray(k%) + Dist * Ray(k% + 3)
        xt(k%) = xh(k%)
    Next k%

    xn(1) = 0: xn(2) = 0: xn(3) = -1
    FlagPlot% = 1

End Sub

Sub Ray.Pixel (Nx%, Ny%, Ray())
    ' Calculates the ray from camera through a given pixel
    AliasX = AntiAlias * (Rnd * 1 - .5)
    AliasY = AntiAlias * (Rnd * 1 - .5)
    For k% = 1 To 3
        Ray(k%) = xCam(k%)
        Ray(k% + 3) = xCenter(k%) - xCam(k%) + ((Nx% + AliasX - xScrCenter%) * UScreen(k%) + (Ny% + AliasY - yScrCenter%) * VScreen(k%)) / ScrHeight
    Next k%
    Ray(7) = Module(Ray(4), Ray(5), Ray(6))
    For k% = 4 To 6: Ray(k%) = Ray(k%) / Ray(7): Next k%
End Sub

Sub Ray.Reflect (xt!(), xn!(), Ray())
    ' Calculate the refelcted ray from a hit point
    Dim xi!(3), Bias!(3)
    For k% = 1 To 3: xi!(k%) = xt!(k%) - xCam(k%): Next k%

    Coeff = 2 * (xi!(1) * xn!(1) + xi!(2) * xn!(2) + xi!(3) * xn!(3))

    For k% = 1 To 3
        Bias!(k%) = 1 + AntiAlias * (Rnd - .5) * .01
        Ray(k%) = xt!(k%)
        Ray(k% + 3) = (xi!(k%) - Coeff * xn!(k%)) * Bias!(k%)
    Next k%

    Ray(7) = Module(Ray(4), Ray(5), Ray(6))
    For k% = 4 To 6: Ray(k%) = Ray(k%) / Ray(7): Next k%
    Ray(7) = .001
End Sub

Sub Ray.Shadow (xt!(), Ray())
    ' Calculate the ray from a hit point to the point of light
    Dim Bias!(3)
    For k% = 1 To 3
        Bias!(k%) = 1 + AntiAlias * (Rnd - .5) * .01
        Ray(k%) = xt(k%) 'xLight(k%)
        Ray(k% + 3) = -(xt!(k%) - xLight(k%)) * Bias!(k%)
    Next k%
    Ray(7) = Module(Ray(4), Ray(5), Ray(6))
    For k% = 4 To 6: Ray(k%) = Ray(k%) / Ray(7): Next k%
    Ray(7) = .01
End Sub

Function Texture.Object! (Id%, FlagPlot%, xh(), xn())
    ' Textures and Bummapping applied to the shapes

    ' Wood : v = 20*P3D, W= v - INt(v)
    ' Marble = M = Cos(x+P3D)

    '==============================
    ' Default Texture: Plain colour
    If shape(Id%).Texture = 0 Then Texture.Object! = shape(Id%).HueRef + 90 * (FlagPlot% - 1)

    '==============================
    ' 2D textures
    If shape(Id%).Texture < 20 Then

        If shape(Id%).Degree = 1 Then
            Select Case shape(Id%).NumType
                Case 100, 101, 104, 201, 204
                    p = 2 * (xh(1) - Int(xh(1))) - 1
                    q = 2 * (xh(2) - Int(xh(2))) - 1
                Case 102, 105, 202, 205
                    p = 2 * (xh(2) - Int(xh(2))) - 1
                    q = 2 * (xh(3) - Int(xh(3))) - 1
                Case 103, 106, 203, 206
                    p = 2 * (xh(1) - Int(xh(1))) - 1
                    q = 2 * (xh(3) - Int(xh(3))) - 1
                Case Else
            End Select
        End If

        If shape(Id%).Degree = 2 Then
            p = Arg(xh(1), xh(2)) / Pi
            q = Arg(Sqr(xh(1) ^ 2 + xh(2) ^ 2), xh(3)) / Pi
        End If

        If shape(Id%).Degree = 3 Then
            p! = xh(1) / shape(Id%).dRef2 / 2
            q! = xh(2) / shape(Id%).dRef2
        End If

        If shape(Id%).Degree = 4 Then
            p = -Arg(Sqr(xh(1) ^ 2 + xh(2) ^ 2), xh(3)) / Pi
            q = Arg(xh(1), xh(2)) / Pi
        End If

        Select Case shape(Id%).Texture
            Case 11
                Texture.Object! = shape(Id%).HueRef + Mandelbrot(p, q)
            Case 12
                Texture.Object! = shape(Id%).HueRef + Julia(p, q)
            Case 13
                Texture.Object! = 720 * Perlin(p, q)
            Case 14
                Noise = Cos(p * 10 + q * 7 + 6 * Perlin!(p, q))
                Texture.Object! = 200 + 320 * Noise
            Case 15
                Texture.Object! = 30 * Cos(10 * p + 10 * Perlin!(p, q))
            Case 16 ' Mandelbrot Ring
                nMdb% = 8
                q = q * nMdb% / 2
                q = (q - Int(q + .5)) / 4
                Texture.Object! = Mandelbrot(2 * p, 2 * q)
            Case 17 ' ChessBoard
                Flagx% = Abs(Int(xh(1) / shape(Id%).dRef1) Mod (2))
                Flagy% = Abs(Int(xh(2) / shape(Id%).dRef2) Mod (2))
                If Flagx% = Flagy% Then
                    Texture.Object! = shape(Id%).HueRef + 90
                Else
                    Texture.Object! = shape(Id%).HueRef
                End If

            Case Else
                Texture.Object! = shape(Id%).HueRef + 90 * (FlagPlot% - 1)

        End Select

    End If

    '========================================
    ' Perlin 3D Hue based texture. For non flat shapes
    If shape(Id%).Texture = 20 Then
        p! = xn(1)
        q! = xn(2)
        r! = xn(3)
        Texture.Object! = 720 * Cos(Pi * (p! + q! + Perlin3D!(p!, q!, r!)))
    End If

    ' ========================================
    ' Sea texture (designed for an horizontal plan)
    If shape(Id%).Texture = 30 Then
        e = 1 / 90 * (1 + .25 * Rnd)
        p! = 2 * (xh(1) / 7 - Int(xh(1) / 7)) - 1
        q! = 2 * (xh(2) / 7 - Int(xh(2) / 7)) - 1
        dx! = Perlin!(p! - e, q!) - Perlin!(p! + e, q!)
        dy! = Perlin!(p!, q! - e) - Perlin!(p!, q! + e)
        xn(1) = 15 * dx! * Cos((xh(1) + xh(2)) * Pi / 2) ^ 2
        xn(2) = 15 * dy!
        xn(3) = 1
        Norm = Module(xn(1), xn(2), xn(3))
        For k% = 1 To 3: xn(k%) = xn(k%) / Norm: Next k%
        Texture.Object! = shape(Id%).HueRef
    End If

    ' ========================================
    ' Sky texture (designed for an horizontal plan)
    If shape(Id%).Texture = 31 Then
        p = 2 * (xh(1) / 100 - Int(xh(1) / 100)) - 1
        q = 2 * (xh(2) / 100 - Int(xh(2) / 100)) - 1
        xn(1) = 5 * Perlin!(p, q): xn(2) = 5 * Perlin!(q, p): xn(3) = -1
        Norm = Module(xn(1), xn(2), xn(3))
        For k% = 1 To 3: xn(k%) = xn(k%) / Norm: Next k%
        Texture.Object! = shape(Id%).HueRef
    End If

    ' ========================================
    ' Nephroid Texture
    If shape(Id%).Texture = 32 Then
        AlphaRef = 55 * Pi / 180
        xp0! = xh(1) * Cos(AlphaRef) + xh(2) * Sin(AlphaRef)
        yp0! = -xh(1) * Sin(AlphaRef) + xh(2) * Cos(AlphaRef)
        If xp0! < 0 Then
            Temp! = Nephroid(xp0!, yp0!, shape(Id%).dRef1)
            If Temp! > 0 Then
                shape(Id%).Brilliance = 1 / (1 + Temp! * 100)
            Else
                shape(Id%).Brilliance = 1 / (1 - Temp! * 1000)
            End If
        Else
            Temp! = Nephroid(0, yp0!, shape(Id%).dRef1)
            shape(Id%).Brilliance = 1 / (1 - Temp! * 1000)
        End If
    End If

    ' ========================================
    ' Bumpmapping (based on Perlin Noises)
    ' ========================================
    ' Bumpmapping for non flat shapes
    If shape(Id%).Bump = 1 Then
        e = .01
        p! = xn(1) * .9
        q! = xn(2) * .9
        r! = xn(3) * .9
        dx! = Perlin3D!(p! - e, q!, r!) - Perlin3D!(p! + e, q!, r!)
        dy! = Perlin3D!(p!, q! - e, r!) - Perlin3D!(p!, q! + e, r!)
        dZ! = Perlin3D!(p!, q!, r! - e) - Perlin3D!(p!, q!, r! + e)
        xn(1) = xn(1) + 5 * dx!
        xn(2) = xn(2) + 5 * dy!
        xn(3) = xn(3) + 5 * dZ!
        Norm = Module(xn(1), xn(2), xn(3))
        For k% = 1 To 3: xn(k%) = xn(k%) / Norm: Next k%
        Texture.Object! = shape(Id%).HueRef + 90 * (FlagPlot% - 1)
    End If

    ' Bumpmapping for flat shapes
    If shape(Id%).Bump = 2 Then

        Select Case shape(Id%).NumType
            Case 101, 104, 201, 204
                p! = xh(1) / shape(Id%).dRef1
                q! = xh(2) / shape(Id%).dRef2
                e = .01
                dx! = Perlin!(p! - e, q!) - Perlin!(p! + e, q!)
                dy! = Perlin!(p!, q! - e) - Perlin!(p!, q! + e)
                xn(1) = xn(1) + 5 * dx!
                xn(2) = xn(2) + 5 * dy!
                xn(3) = xn(3)

            Case 102, 105, 202, 205
                p! = xh(2) / shape(Id%).dRef1
                q! = xh(3) / shape(Id%).dRef2
                e = .01
                dx! = Perlin!(p! - e, q!) - Perlin!(p! + e, q!)
                dy! = Perlin!(p!, q! - e) - Perlin!(p!, q! + e)

                xn(2) = xn(2) + 5 * dx!
                xn(3) = xn(3) + 5 * dy!
                xn(1) = xn(1)

            Case 103, 106, 203, 206
                p! = xh(1) / shape(Id%).dRef1
                q! = xh(3) / shape(Id%).dRef2
                e = .01
                dx! = Perlin!(p! - e, q!) - Perlin!(p! + e, q!)
                dy! = Perlin!(p!, q! - e) - Perlin!(p!, q! + e)
                xn(1) = xn(1) + 5 * dx!
                xn(3) = xn(3) + 5 * dy!
                xn(2) = xn(2)
            Case Else
        End Select

        Norm = Module(xn(1), xn(2), xn(3))
        For k% = 1 To 3: xn(k%) = xn(k%) / Norm: Next k%
    End If

    If shape(Id%).Bump = 3 Then
        ' Water effects
        Amplitude = .1
        Fade = .2
        Omega = 8
        r = Sqr(xh(1) ^ 2 + xh(2) ^ 2)
        'A = 1 / (1 + Fade * R) * SIN(Omega * R)

        dAdR = -Fade / (1 + Fade * r) * Sin(Omega * r) + Omega / (1 + Fade * r) * Cos(Omega * r)

        If r <> 0 Then
            dRdX = xh(1) / r
            drdY = xh(2) / r
        Else
            dRdX = 0: drdY = 0
        End If

        dadX = Amplitude * dAdR * dRdX
        dadY = Amplitude * dAdR * drdY
        xn(1) = -dadX: xn(2) = -dadY: xn(3) = 1
        Norm = Module(xn(1), xn(2), xn(3))
        For k% = 1 To 3: xn(k%) = xn(k%) / Norm: Next k%
    End If

End Function