S.15　Double Toruses (Tori)

Public Sub Proc15() 'Double torus
    Dim id
    Dim r1, r2, k1, k2, ip, lastp, zzzw, ippw
    Dim th1, th2, th3, xs1, ys1, xs2, ys2, xs3, ys3, xs4, ys4
    Dim thz, thy, x1, y1, zw, xw
    Dim xs(50, 17), ys(50, 17), zs(50, 17), zc(50, 17)
    Dim zzz(), ipp()
    Dim points(4) As POINTAPI
    r1 = 150: r2 = 50: k1 = 24: k2 = 16
    lastp = 2 * k1 * k2
    ReDim zzz(lastp), ipp(lastp)
    Select Case Ivar
        Case 1: thy = 0.2: th3 = 0
        Case 2: thy = 0.6: th3 = 0
        Case 3: thy = 0.6: th3 = 1
        Case 4: thy = -0.6: th3 = 1
        Case Else: thy = 0.2: th3 = 0
    End Select
    Dperas
    Dpwind 0, 0, 640
'------ Camera Position
    CamP(1) = 1200
    CamP(2) = 0
    CamP(3) = 0
    Dpcam 0.6, 0, Flen
    Dpwind 0, 0, 640
    N = 0
    For id = -1 To 1 Step 2
        For th1 = 0 To 2 * pie + 0.1 Step 2 * pie / k1
            N = N + 1: m = 0
            For th2 = 0 To 2 * pie + 0.1 Step 2 * pie / k2
                m = m + 1
                Y = r1 + r2 * Cos(th2)
                z = r2 * Sin(th2)
                X = 0
                thz = th1
                Rotz X, Y, X, Y, thz
                Y = Y + r1 / 2 * id
                If id = 1 Then
                    zw = z: xw = X
                    z = zw * Cos(pie / 2) - xw * Sin(pie / 2)
                    X = zw * Sin(pie / 2) + xw * Cos(pie / 2)
                End If
                thz = th3
                Rotz X, Y, X, Y, thz
                
                Roty z, X, z, X, thy
                
                Projection X, Y, z
                ys(N, m) = Gxxx
                zs(N, m) = Gyyy
                xs(N, m) = X
            Next th2
        Next th1
    Next id
'------ distance from the camera
    ip = 1
    For N = 1 To k1 * 2 + 1
        If N <> k1 + 1 Then
            For m = 1 To k2
                zc(N, m) = (xs(N, m) + xs(N + 1, m + 1)) / 2
                zzz(ip) = zc(N, m): ipp(ip) = ip
                ip = ip + 1
            Next m
        End If
    Next N
'sort---------------------------------
    For i = 2 To lastp
        For j = i - 1 To 1 Step -1
            If zzz(j) > zzz(j + 1) Then
                zzzw = zzz(j): zzz(j) = zzz(j + 1): zzz(j + 1) = zzzw
                ippw = ipp(j): ipp(j) = ipp(j + 1): ipp(j + 1) = ippw
            End If
        Next j
    Next i
'draw--------------------------------
    For ip = 1 To lastp
            N = Int((ipp(ip) - 1) / k2) + 1
            If N > k1 Then N = N + 1
            m = (ipp(ip) - 1) Mod k2 + 1
            xs1 = ys(N, m): ys1 = zs(N, m)
            xs2 = ys(N + 1, m): ys2 = zs(N + 1, m)
            xs3 = ys(N + 1, m + 1): ys3 = zs(N + 1, m + 1)
            xs4 = ys(N, m + 1): ys4 = zs(N, m + 1)
            If Int((ys4 - ys2) * (xs3 - xs1) + _
                 (ys3 - ys1) * (xs2 - xs4)) > 100 Then
                points(1).X = Ipixelx(xs1)
                points(2).X = Ipixelx(xs2)
                points(3).X = Ipixelx(xs3)
                points(4).X = Ipixelx(xs4)
                points(1).Y = Ipixely(ys1)
                points(2).Y = Ipixely(ys2)
                points(3).Y = Ipixely(ys3)
                points(4).Y = Ipixely(ys4)
                Dpolyg points, 4
            End If
    Next ip
    Dptext 150, 200, ("thy=" + CStr(thy))
End Sub