open GeometRi

open System

​

printfn "Ptolemy's construction of a pentagon"

printfn "Draw a regular pentagon inside the circle"

printfn ""

​

// We will follow the notation from 

// https://www.cut-the-knot.org/pythagoras/pentagon.shtml

// Read this article for detailed description of the algorithm

​

// Create circle with center in point O and choose arbitrary point C on circle

// We will start our pentagon from this point

let O = Point3d(0.0, 0.0, 0.0)

let circle = Circle3d(O, 1.0, Vector3d(0.0, 0.0, 1.0))

let C = circle.ParametricForm(0.0)

​

// Draw a diameter through point C (point D is the opposite point of the diameter)

let l = Line3d(C, O)

let CD = l.IntersectionWith circle |> unbox<Segment3d>

​

// Need to check if we get the correct point

let D = if CD.P1 = C then CD.P2 else CD.P1

​

// Draw a orthogonal diameter AB

​

// First, draw two circles with centers in points C and D

// The intersection of these circles gives a line, orthogonal to the segment CD

let circle1 = Circle3d(C, CD.Length, circle.Normal)

let circle2 = Circle3d(D, CD.Length, circle.Normal)

let line = (circle1.IntersectionWith circle2 |> unbox<Segment3d>).ToLine

​

// Now find the intersection of line with our circle

// This gives points A and B

let AB = line.IntersectionWith circle |> unbox<Segment3d>

let A = AB.P1

let B = AB.P2

​

// Find the center point of the segment OB (this will be point M)

let M = (O + B) / 2.0

​

// Draw a circle with center in point M and radius MC

let radius = M.DistanceTo C

let tmp_circle0 = Circle3d(M, radius, circle.Normal)

​

// Find the intersection with diameter AB (this will be point N)

let ON = tmp_circle0.IntersectionWith AB |> unbox<Segment3d>

let N = if ON.P1 = O then ON.P2 else ON.P1

​

// CN is the length of the pentagon's side

let side_len = C.DistanceTo N 

​

// The first point of the pentagon is point C

let P1 = C

​

// Draw a circle from point C with radius 'side_len' and this gives us two more points

let tmp_circle1 = Circle3d(C, side_len, circle.Normal)

let tmp_segment = tmp_circle1.IntersectionWith circle |> unbox<Segment3d>

let P2 = tmp_segment.P1

let P5 = tmp_segment.P2

​

// Repeat this two more times to get remaining points

let tmp_circle2 = Circle3d(P2, side_len, circle.Normal)

let tmp_segment2 = tmp_circle2.IntersectionWith circle |> unbox<Segment3d>

let P3 = if tmp_segment2.P1 = P1 then tmp_segment2.P2 else tmp_segment2.P1

​

let tmp_circle3 = Circle3d(P3, side_len, circle.Normal);

let tmp_segment3 = tmp_circle3.IntersectionWith circle |> unbox<Segment3d>

let P4 = if tmp_segment3.P1 = P2 then tmp_segment3.P2 else tmp_segment3.P1 

​

// Points P1, P2, P3, P4 and P5 define our pentagon

// Let's check if interion angles are equal 108 degrees (with accuracy 1e-8)

let angle1 = Vector3d(P1, P2).AngleToDeg(Vector3d(P1, P5))

let angle2 = Vector3d(P2, P1).AngleToDeg(Vector3d(P2, P3))

let angle3 = Vector3d(P3, P2).AngleToDeg(Vector3d(P3, P4))

let angle4 = Vector3d(P4, P3).AngleToDeg(Vector3d(P4, P5))

let angle5 = Vector3d(P5, P4).AngleToDeg(Vector3d(P5, P1))

​

let tolerance = 1e-8;

if Math.Abs(angle1 - 108.0) < tolerance &&

   Math.Abs(angle2 - 108.0) < tolerance &&

   Math.Abs(angle3 - 108.0) < tolerance &&

   Math.Abs(angle4 - 108.0) < tolerance &&

   Math.Abs(angle5 - 108.0) < tolerance then

    printfn "Your construction is correct!"

    printfn "%A" P1

    printfn "%A" P2

    printfn "%A" P3

    printfn "%A" P4

    printfn "%A" P5