#!/usr/bin/perl -CSDA

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

# Calculate the radius of a circle that circumscribes a set of adjacent cords
# described by chord lengths. Such a polygon is called a cyclic polygon.
# The radius is called the circumradius.

# This is being used to do test research for an Openscad project that needs
# this.

# Closed functions to calculate a circumradius given a list of chord
# lengths have been derived in the general case up to heptagons. These
# aren't really practical though.

# The plan is to handle three cases on whether twice the circumradius is
# less, equal or greater than the longest chord, using similar processes
# with some tweaks. For here at least, I want to try using Halley's
# irrational method for iterating on the circumradius.

# Each chord will use an angle of 2*asin(C/(2*R)) radians.
# The total of angles needs to add up to 2*pi().
# The derivative (dR) of 2*arcsine(C/(2*R)) is
# -2*C/R/sqrt(4*R*R-C*C)
# The second derivative (d2R) is
# 2*C/(R*R)/sqrt(4*R*R-C*C)*(2+C*C/(4*R*R-C*C))
# If the center of the circumscribed circle is outside of the cyclic polygon,
# then we effectively use that the arc of the largest side is equal
# to the sum of the other arcs, rather than that the sum is 2 pi.

use utf8;
use utf8::all;
use strict;
use warnings;
use Scalar::Util qw(looks_like_number);
use Math::Trig;

# There should be a least three arguments with some relative size restrictions.
# And all values should be numeric.
#
my @chords = @ARGV;

if (scalar @chords < 3) {
  die "There must be at least three chord lengths.\n";
}

my $sum = 0;
my $max = 0;
my $imax = 0;
my $i = 0;
foreach my $chord (@chords) {
  if (!looks_like_number($chord)) {
    die "The expected arguments are numbers.\n";
  }
  if ($chord <= 0) {
    die "The expected arguments are positive numbers.\n";
  }
  $sum += $chord;
  if ($max < $chord) {
    $max = $chord;
    $imax = $i;
  }
  $i++;
}

foreach my $chord (@chords) {
  if (2*$chord >= $sum) {
    die "One of the chords is too long.\n";
  }
}

# See which case we are in.
# We do this by assuming that the center of the circumscribed circle is
# inside the cyclic polygon and then test if the radius is shorter than
# half he longest chord. If it is smaller, then there is a contradicton
# and the center must be outside of the cyclic polygon.

my $target = 2*pi();
my $r = $max/2;
my $arc = 0;
foreach my $chord (@chords) {
  $arc += 2*asin($chord/(2*$r));
}
if ($arc == $target) {
  print "Circumradius = ", $r, "\n";
  exit;
} elsif ($arc < $target) {
  $target = 0;
  $chords[$imax] = -$chords[$imax];
}

# The derivative of asin is ill conditioned when the radius is half
# the chord length, so we need to move away from that, before
# starting to iterate. It will always be bigger, so move in that
# direction. We also need to stay near half the longest chord or we
# may get outside of an area where we will converge.

$r *= 1.001;

for (my $i = 0; $i <= 10; $i++) {
  my $arc = 0;
  my $der = 0;
  my $der2 = 0;
  foreach my $chord (@chords) {
    $arc += 2*asin($chord/(2*$r));
    $der += -2*$chord/sqrt(4*$r*$r-$chord*$chord)/$r;
    $der2 += 2*$chord/($r*$r)/sqrt(4*$r*$r-$chord*$chord)*(2+$chord*$chord/(4*$r*$r-$chord*$chord));
  }
  print "Circumradius = ", $r, "  Angle = ", $arc, "\n";
  $r = $r+2*($target-$arc)/$der/(1+sqrt(1-2*($target-$arc)*$der2/($der*$der)));
}