Imports GeometRi

Imports System

​

Public Module Module1

​

    Public Sub Main()

​

        Console.WriteLine("Rytz's construction example")

​

        ' We will follow the notation from 

        ' https://en.wikipedia.org/wiki/Rytz%27s_construction

        ' Read this article for detailed description of the algorithm

​

        ' Define arbitrary circle And plane

        Dim circle As Circle3d = New Circle3d(New Point3d(2, 3, -1), 5, New Vector3d(1, 1, 1))

        Dim plane As Plane3d = New Plane3d(New Point3d(-1, -5, 0), New Vector3d(3, 0, 1))

​

        ' Get two points on the circle at orthogonal diameters

        Dim d1 As Point3d = circle.ParametricForm(0)

        Dim d2 As Point3d = circle.ParametricForm(Math.PI / 2)

​

        ' Project onto plane

        Dim C As Point3d = circle.Center.ProjectionTo(plane)

        Dim P As Point3d = d1.ProjectionTo(plane)

        Dim Q As Point3d = d2.ProjectionTo(plane)

​

        ' Rotate point P by 90 degrees

        Dim r As Rotation = New Rotation(plane.Normal, Math.PI / 2)

        Dim P_prime As Point3d = P.Rotate(r, C)

​

        ' Determine the center D of the line segment P'-Q

        Dim D As Point3d = (P_prime + Q) / 2

​

        ' Draw the circle with center D throuugh point C

        Dim radius As Double = D.DistanceTo(C)

        Dim circle_DC As Circle3d = New Circle3d(D, radius, plane.Normal)

​

        ' Intersect the circle And the line P'-Q. Intersection points are A and B.

        ' (some additional treatment is needed when circle is parallel to plane,

        '  not going to do it here)

        Dim line As Line3d = New Line3d(P_prime, Q)

        Dim AB As Segment3d = CType(line.IntersectionWith(circle_DC), Segment3d)

        Dim A As Point3d = AB.P1

        Dim B As Point3d = AB.P2

​

        ' Find the length of semi-axes

        Dim semi_axis_a As Double = Q.DistanceTo(A)

        Dim semi_axis_b As Double = Q.DistanceTo(B)

​

        ' Find the direction of semi-axes

        Dim vec_CA As Vector3d = New Vector3d(C, A)

        Dim vec_CB As Vector3d = New Vector3d(C, B)

        ' renormalize vectors

        vec_CA = semi_axis_b * vec_CA.Normalized

        vec_CB = semi_axis_a * vec_CB.Normalized

​

        ' Now, construct the projection

        Dim My_Ellipse As Ellipse = New Ellipse(C, vec_CA, vec_CB)

​

        ' Compare with more efficient numerical algorithm

        If (My_Ellipse = circle.ProjectionTo(plane)) Then

            Console.WriteLine("Your construction is correct!")

            Console.WriteLine(My_Ellipse.ToString())

        End If

​

​

        Console.ReadLine()

    End Sub

​

End Module