Sub EnumTetraOpt()
    Dim j As Integer, l(1 To 6) As Integer
    Dim n1 As Integer, n2 As Integer, n3 As Integer, n4 As Integer, n5 As Integer, n6 As Integer
    Dim k As Integer, ix As Integer, iy As Integer
    Dim cb As Integer, pm As Integer
    Dim n As Long, nT As Long, nD As Long, nI As Long
    Dim nMax As Long, nPrint As Long, nTot As Double
    Dim Comb63(1 To 20) As String * 3
    Dim Prm3(1 To 6) As String * 3
    Dim i(0 To 100) As Integer

    Dim nComb As Double
    Dim comb As Double, temp As Double
    
    Comb63(1) = "123"
    Comb63(2) = "124"
    Comb63(3) = "125"
    Comb63(4) = "126"
    Comb63(5) = "134"
    Comb63(6) = "135"
    Comb63(7) = "136"
    Comb63(8) = "145"
    Comb63(9) = "146"
    Comb63(10) = "156"
    Comb63(11) = "234"
    Comb63(12) = "235"
    Comb63(13) = "236"
    Comb63(14) = "245"
    Comb63(15) = "246"
    Comb63(16) = "256"
    Comb63(17) = "345"
    Comb63(18) = "346"
    Comb63(19) = "356"
    Comb63(20) = "456"
    'Comb63(k) complementary to Comb63(21-k)
    
    Prm3(1) = "123"
    Prm3(2) = "132"
    Prm3(3) = "213"
    Prm3(4) = "231"
    Prm3(5) = "312"
    Prm3(6) = "321"

On Error GoTo handlerr

    Debug.Print Time$
    For nMax = 1 To 200
        i(0) = 0
        nTot = 0: nT = 0
        For k = 1 To nMax: i(k) = 6: Next k
            
        Do
            comb = Fact(6)
            For k = 1 To nMax:
                temp = i(k) - i(k - 1)
                comb = comb / Fact(temp)
            Next k
            nTot = nTot + comb
            nComb = nComb + 1
            ix = 1
            Do
                DoEvents
                j = 1: iy = 1
                For k = 1 To nMax + 1
                    For j = 1 To i(k) - i(k - 1)
                        l(iy) = k
                        iy = iy + 1
                    Next j
                Next k
                
                For cb = 1 To 20
                    n1 = Val(Mid$(Comb63(cb), 1, 1))
                    n2 = Val(Mid$(Comb63(cb), 2, 1))
                    n3 = Val(Mid$(Comb63(cb), 3, 1))
                    For pm = 1 To 6
                        n4 = Val(Mid$(Comb63(21 - cb), Val(Mid$(Prm3(pm), 1, 1)), 1))
                        n5 = Val(Mid$(Comb63(21 - cb), Val(Mid$(Prm3(pm), 2, 1)), 1))
                        n6 = Val(Mid$(Comb63(21 - cb), Val(Mid$(Prm3(pm), 3, 1)), 1))
                        
                        If Triangular(l(n1), l(n2), l(n3)) And Triangular(l(n1), l(n4), l(n5)) _
                            And Triangular(l(n2), l(n5), l(n6)) And Triangular(l(n3), l(n6), l(n4)) Then
                            If (l(n6) < dmax(l(n1), l(n2), l(n3), l(n5), l(n4))) Then
                                If l(n6) > dmin(l(n1), l(n2), l(n3), l(n5), l(n4)) Then
                                    nT = nT + comb
                                    GoTo nxt
                                End If
                            End If
                        End If
                    Next pm
                Next cb
nxt:
                
                If i(ix) > 0 Then
                    i(ix) = i(ix) - 1
                    ix = 1
                    Exit Do
                Else
                    Do
                        ix = ix + 1
                        If ix = nMax Then Exit Do
                    Loop Until i(ix) > 0
                    If ix = nMax Then Exit Do
                    i(ix) = i(ix) - 1
                    For k = 1 To ix - 1: i(k) = i(ix): Next k
                    ix = 1
                    Exit Do
                End If
            Loop
        Loop While ix < nMax
        Debug.Print Time$, "n="; nMax, "nTet="; nT
    Next nMax
    Beep
    End
    
handlerr: Stop: Resume Next
End Sub

Function Triangular(x As Integer, y As Integer, z As Integer) As Boolean
    Triangular = Sgn(2 * (x ^ 2 * y ^ 2 + y ^ 2 * z ^ 2 + z ^ 2 * x ^ 2) - (x ^ 4 + y ^ 4 + z ^ 4)) = 1
End Function

Function Fact(x) As Double
Dim k As Integer
    Fact = 1
    For k = 2 To x
        Fact = Fact * k
    Next k
End Function

Function dmax(a As Integer, b As Integer, c As Integer, d As Integer, e As Integer) As Double
    'triangles abc, ade with a in common on the same side
    'evaluates distance of vertices bc,de

    Dim d2 As Double
    
    d2 = ((a ^ 2 + b ^ 2 - c ^ 2) / (2 * a) - (a ^ 2 + d ^ 2 - e ^ 2) / (2 * a)) ^ 2 + _
        ((2 * a ^ 2 * b ^ 2 + 2 * a ^ 2 * c ^ 2 + 2 * b ^ 2 * c ^ 2 - a ^ 4 - b ^ 4 - c ^ 4) ^ (1 / 2) / (2 * a) + _
          (2 * a ^ 2 * d ^ 2 + 2 * a ^ 2 * e ^ 2 + 2 * d ^ 2 * e ^ 2 - a ^ 4 - d ^ 4 - e ^ 4) ^ (1 / 2) / (2 * a)) ^ 2
    dmax = d2 ^ (1 / 2)
End Function

Function dmin(a As Integer, b As Integer, c As Integer, d As Integer, e As Integer) As Double
    'triangles abc, ade with a in common on opposite sides
    'evaluates distance of vertices bc,de

    Dim d2 As Double
    
    d2 = ((a ^ 2 + b ^ 2 - c ^ 2) / (2 * a) - (a ^ 2 + d ^ 2 - e ^ 2) / (2 * a)) ^ 2 + _
        ((2 * a ^ 2 * b ^ 2 + 2 * a ^ 2 * c ^ 2 + 2 * b ^ 2 * c ^ 2 - a ^ 4 - b ^ 4 - c ^ 4) ^ (1 / 2) / (2 * a) - _
          (2 * a ^ 2 * d ^ 2 + 2 * a ^ 2 * e ^ 2 + 2 * d ^ 2 * e ^ 2 - a ^ 4 - d ^ 4 - e ^ 4) ^ (1 / 2) / (2 * a)) ^ 2
    dmin = d2 ^ (1 / 2)
End Function