// SPDX-License-Identifier: LicenseRef-PD-hp OR CC0-1.0 OR 0BSD OR MIT-0 OR MIT

// Return the circumradius for a list of chord lengths and a list of
// angles. When the center of the circle is outside the cyclic polygon,
// the angle for the longest side needs to be complimemted to save
// special handling later.

// There isn't a general closed form solution for this. So we look at
// how the total arc length changes with respect to he radius and use
// Halley's irrational method to iterate to a result.

// The approach needs to be tweaked depending on whether the center
// of the circle is inside, outside or on the border of the cyclic
// polygon. (As above, this isn't affected by the order of the chords.)

function circumradius(chord_lengths) = 

  assert(is_list(chord_lengths), "A list of chord lengths was not provided.")
  assert(len(chord_lengths)>=3, "There must be at least three chords.")
// is_num_v is defined in utils.scad.
  assert(is_num_v(chord_lengths), "Not all lengths are numbers.")
  assert(min(chord_lengths)>0, "Not all lengths are positive.")
// sum is defined in utils.scad.
  assert(max(chord_lengths)<(sum(chord_lengths)/2), "A polygon can't be made from the chords.")

// Trying the circumradius equal to half the longest chord is the critical
// point for whether or not the center of the circle is inside or outside
// the polygon.
  let (r=max(chord_lengths)/2)
  let (arc=2*sum([for (chord_length=chord_lengths) asin(chord_length/(2*r))]))

// We need to start near, but not at a radius equal to half the longest chord.
  arc>360?
    iterate(r*1.001, chord_lengths, 360):
    arc<360?
// For the outside the polygon case, there can be only one max length chord.
    iterate(r*1.001, [for (chord_length=chord_lengths) chord_length==max(chord_lengths)?-chord_length:chord_length], 0):
    [r, [for (chord_length=chord_lengths) 2*asin(chord_length/(2*r))]]
;

// Use Halley's irrational method to converge on the answer.
// Target will be 360 for the inside case or 0 for the outside case.

function iterate(r, chord_lengths, target) =

  let (arc=2*sum([for (chord_length=chord_lengths) asin(chord_length/(2*r))]))
  let (der=-360/PI*sum([for (chord_length=chord_lengths) chord_length/sqrt(4*r*r-chord_length*chord_length)/r]))
  let(der2=360/PI*sum([for (chord_length=chord_lengths) chord_length/(r*r)/sqrt(4*r*r-chord_length*chord_length)*(2+chord_length*chord_length/(4*r*r-chord_length*chord_length))]))
  let(newr=r+2*(target-arc)/der/(1+sqrt(1-2*(target-arc)*der2/(der*der))))
  abs(newr-r)/r>1e-9?
    iterate(newr, chord_lengths, target):
    [newr, [for (chord_length=chord_lengths) chord_length<0?360+2*asin(chord_length/(2*newr)):2*asin(chord_length/(2*newr))]]
;