#!/bin/perl -w
use strict;

my($squares, $Klein, $narrow, $dump, $more, $guidelines) = 5;
$Klein = shift if ($ARGV[0] || 0) eq 'Klein';
$narrow = shift  if ($ARGV[0] || 0) eq 'narrow';
$dump = shift  if ($ARGV[0] || 0) eq 'dump';
$guidelines = shift  if ($ARGV[0] || 0) eq 'guidelines';
$more = shift  if ($ARGV[0] || 0) eq 'more';
$squares = shift  if @ARGV;
die "Usage: 
  $0 [Klein] [narrow] [dump] [guidelines] [more] [N]
     N defaults to 5\n" if @ARGV;

=pod

Calculating a square grid on Poincaré circle with angles 72°: diagonals break
the square into triangles with angles 36°, 36° and 90°.  Putting a sharp angle
at the center O, denote by A the other sharp angle.  The other side at A is a 
Euclidean circle c with center O', and ∠OAO' is 36°+90°.  Hence (Euclidean
distances) |AO'|/sin 36° = |OA|/sin 18°, |OO'|/cos 36° = |OA|/sin 18°.  Since
the length of tangent to c from O is 1, |OO'|²-|AO'|²=1.  Hence
  |OA|²/sin² 18° · (cos² 36° - sin² 36°) = 1,
or |OA|² = sin² 18°/(cos² 36° - sin² 36°) = sin² 18°/cos 72° = sin² 18°/sin 18°,
and |OA| = √sin 18°.

More generally, for angle α (instead of 36°)
  |AO'|/sin α = |OA|/cos 2α, |OO'|/cos α = |OA|/cos 2α,
  |OA|²/cos² 2α · (cos² α - sin² α) = 1,
  |OA| = √cos 2α

Conformal transformations of a circle
=====================================

To send z₀ to 0:  w = (z-z₀)/(1 - z̅₀ z) (preserved diameter through z₀), or
 w = -1/z̅₀ (1 + (1/z̅₀ - z₀)/(z-1/z̅₀))
 -(1/z̅₀ - z₀)/(1 + z̅₀w) + 1/z̅₀ = z
 z = (w+z₀)/(1+z̅₀w).

=cut

use Math::Pari qw(PARI PARImat PARIcol I Pi real imag conj setprecision);
setprecision(120);

sub remove_nl ($) { (my $in = shift) =~ s/\n//g; $in }
sub parse_as_gp ($) {
  my $in = shift;
  $in =~ s/(\\[^\\]|[^"\s\\]|\n|"([^"\\]|\\.)*")|[^\S\n\r]+|\\\\.*/ defined($1) ? $1 : '' /ges;
  # Now all unneeded whitespace (except LF) and comments are removed
  $in =~ s/^\{(.*?)^}$/ remove_nl $1 /gems;
  $in =~ s/=\n/=/gsm;
  $in =~ s/\\\n//g;
  warn "in: <<$in>>\n";
  my @in = split /\n/, $in;
  $in = pop @in;
  PARI $_ for @in;	# Void context (ignored???)
  PARI $in
}

my $ignore = parse_as_gp <<EOP;
send_to0(z0,z) =
  (z-z0)/(1 - conj(z0)*z)
send_0to(z0,w) =
  (w+z0)/(1 + conj(z0)*w)
relative_rotate(z0,z,factors) =
{  local(zz);
   zz = send_to0(z0,z);
   vector(length(factors), K, send_0to(z0,zz*factors[K]))
}
Factors = vector($squares-1,K,exp(2*Pi*I/$squares*K))
more_points(z0,z) =
  relative_rotate(z0, z, Factors)
EOP

############################

#my $seed = sqrt cos(2*Pi/$squares);			# worse precision???  0.0 vs -0.0
my $seed = sqrt sin(Pi*((1 - PARI(4)/$squares)/2));
sub seed() { $seed }			# make visible from PARI
my @level1 = ($seed, @{ PARI 'more_points(0,seed())' });
my @levels = ([PARI 0], \@level1);

my $hash_pattern = $more ? '%.3f %.3f' : '%.2f %.2f';
sub hash_point($) {my $in = shift; sprintf $hash_pattern, real $in, imag $in}
sub hash_edge($$) {my($p,$pp) = (shift, shift); "$p, $pp"}

sub eps() {1e-14}
sub epseps() {1e-100}
sub find_center($$) {				# find an arc perpendicular to |z|=1 through points
  my($p1,$p2) = (shift, shift);
  return if eps > abs imag($p1*conj $p2);	# not an arc
#  my($pp1,$pp2) = map 1/conj $_, $p1, $p2;	# inversion images, also on arc
#  my $M = 1/PARImat([imag $p1, imag $p2], [-real $p1, -real $p2]);
  my $M = 1/PARImat([real $p1, real $p2], [imag $p1, imag $p2]);
  my($v1,$v2) = map 0.5*(1 + abs($_)**2), $p1, $p2;	# pairings of perp(z0) with z0, 1/bar(z0) are abs(z0)**2, 1
  PARI(1,I) * $M * PARIcol $v1, $v2;
}
sub make_arc($$) {
  my ($p1,$p2) = (shift, shift);
  if ($Klein or not (my($c) = find_center $p1, $p2)) {
    if ($Klein) { 
      $_ = $_ * 2/(1+abs($_)**2) for $p1, $p2;
    }
    sprintf "%f %f moveto %f %f lineto\n", real $p1, imag $p1, real $p2, imag $p2;
  } else {
    my $r = abs($c - $p1);
    # draw counterclockwise; reorder points if needed
    ($p1, $p2) = ($p2, $p1) if 0 > imag(($p2-$c) * conj($p1-$c));
    my($a1, $a2) = map 180/Pi * imag log($_ - $c), $p1, $p2;
    sprintf "%f %f moveto %f %f %f %f %f arc\n", real $p1, imag $p1, real $c, imag $c, $r, $a1, $a2;
  }
}

my %points = map +(hash_point $_, $_), 0, @level1;
my %edges = map +(hash_edge(0, $_), 1), @level1;
my @edges = map [0, $_], @level1;
my $just_reached = [ map [$_, 0], @level1 ];
my(@level2, $to, $from, @onguide) = map real $_, @level1[0..int(($squares+1)/2)];
sub from()	{$from}
sub to()	{$to}
for my $r (2..14) {
  my @reached;
  for my $R (@$just_reached) {
    ($to, $from) = @$R;
    for my $p ( @{ PARI 'more_points(to(),from())' } ) {
      push @edges, [$to, $p];
      warn "unexpected repeated edge: [$to, $p]" if $edges{hash_edge $to, $p} or $edges{hash_edge $p, $to};
      $edges{hash_edge $to, $p}++;
      my $h = hash_point $p;
      # @level2 may contain duplicates???
      unless (exists $points{$h}) {
        imag $p > - eps and push @level2, real $p if $r == 2;
        $points{$h} = $p;
        push @reached, [$p, $to];
        if ($guidelines) {
          for my $X (@level2) {
            push(@onguide, $h), last if epseps > abs($X - real $p);
          }
        }
      }
    }
  }
  warn "Saturated: $r" and last unless @reached;
  $just_reached = \@reached;
}

# Recalculate from Poincare model to Klein one: tan x |---> tan 2x (hyperbolic version)
#@edges = map [map $_ * 2/(1+abs($_)**2), @$_], @edges if $Klein;

if ($dump) {
  print "edge: [@$_]\n" for @edges;
  exit;
}

my $w = $narrow ? 0.0001 : 0.0005;
print <<EOP;
%!
250 250 scale
1.2 1.2 translate
$w setlinewidth
EOP

if ($guidelines) {
  print <<EOP;
0 0 1 setrgbcolor

EOP
  printf <<EOP, real $_, imag $_, $w*3 for map $points{$_}, @onguide;
%f %f %f 0 360 arc closepath stroke
EOP
  print <<EOP;
0 setgray

EOP
}

print make_arc $_->[0], $_->[1] for @edges;

if ($guidelines) {
  print <<EOP;
stroke
1 0 0 setrgbcolor

EOP
  printf <<EOP, $_, $_ for @level2;
%f -1 moveto
%f  1 lineto
EOP
}

print <<EOP;
stroke
showpage
EOP