using GeometRi;

using System;

​

​

public class Program

{

    public static void 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

        Circle3d circle = new Circle3d(new Point3d(2, 3, -1), 5, new Vector3d(1, 1, 1));

        Plane3d plane = new Plane3d(new Point3d(-1, -5, 0), new Vector3d(3, 0, 1));

​

​

        // Get two points on the circle at orthogonal diameters

        Point3d d1 = circle.ParametricForm(0);

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

​

        //Project onto plane

        Point3d C = circle.Center.ProjectionTo(plane);

        Point3d P = d1.ProjectionTo(plane);

        Point3d Q = d2.ProjectionTo(plane);

​

        // Rotate point P by 90 degrees

        Rotation r = new Rotation(plane.Normal, Math.PI / 2);

        Point3d P_prime = P.Rotate(r, C);

​

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

        Point3d D = (P_prime + Q) / 2;

​

        // Draw the circle with center D throuugh point C

        double radius = D.DistanceTo(C);

        Circle3d circle_DC = 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)

        Line3d line = new Line3d(P_prime, Q);

        Segment3d AB = (Segment3d)line.IntersectionWith(circle_DC);

        Point3d A = AB.P1;

        Point3d B = AB.P2;

​

        // Find the length of semi-axes

        double semi_axis_a = Q.DistanceTo(A);

        double semi_axis_b = Q.DistanceTo(B);

​

        // Find the direction of semi-axes

        Vector3d vec_CA = new Vector3d(C, A);

        Vector3d vec_CB = 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

        Ellipse ellipse = new Ellipse(C, vec_CA, vec_CB);

​

        // Compare with more efficient numerical algorithm

        if (ellipse == circle.ProjectionTo(plane))

        {

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

            Console.WriteLine(ellipse.ToString());

        }

​

        Console.ReadLine();

​

​

    }

}

​

​

​