### Poincaré disk model #1

This is my second attempt at drawing the Poincaré disk model in a turtle.

Variation: Poincaré disk model #1 (variation)

Useful Resources:
en.wikipedia.org/wiki/poincar%c3%a9_disk_model
strauss.hosted.uark.edu/papers/hypcomp.pdf
mathcs.clarku.edu/~djoyce/poincare/

#Poincare

```// Poincaré disk model #1. Created by Reinder Nijhoff 2021 - @reindernijhoff
//
// https://turtletoy.net/turtle/d176924430
//
// This is my second attempt at drawing the Poincaré disk model in a turtle. In my first attempt, I tried to follow the
// construction as explained in [1], but failed at deeper recursion.
// Now I am using the same algorithm as David E. Joyce[2]. First the center polygon is constructed, then adjacent polygons
// are constructed by (recursive) reflection.
//
// Geodesics are represented as vec3's u, v, r2 (following  @matt_zucker[3]):
//   * either a circle centered at (u, v) with squared radius r if r2 > 0
//   * or a diameter of the unit circle with normal (u, v) if r2 == 0
//
// [1] https://strauss.hosted.uark.edu/papers/hypcomp.pdf
// [2] https://mathcs.clarku.edu/~djoyce/poincare/
//

Canvas.setpenopacity(.5);

const turtle = new Turtle();
const n = 3; // min=3, max=13, step=1
let   k = 3; // min=3, max=13, step=1
const recursion=7; // min=1, max=10, step=1
const size = 95; // min=10, max=100, step=1
const centerOffset = 0; // min=0, max=1, step=0.01
const drawFullGeodesics = 0; // min=0, max=1, step=1 (No, Yes)

if ( n==3 )      { k = Math.max(k,7); }
else if (n == 4) { k = Math.max(k,5); }
else if (n < 7)  { k = Math.max(k,4); }
else             { k = Math.max(k,3); }

// Hyperbolic geometry functions

// Construction 1.1. Construct a circle through three given non-collinear points A,B, C.
function circleConstruct3Points(a, b, c) {
const lab = bisector(a, b);
const lbc = bisector(b, c);
const center = cross(lab, lbc);
return compass(scale(center,1/center[2]), a);
}

// Construction 1.2. Invert a point through a circle C with center O.
function circleInvert(c, p) {
const po = sub(p, c);
return add(c, scale(po, c[2] / dot(po, po)));
}

// Construction 1.6. Given a circle C with center O, and point A in the exterior of C,
// construct the unique circle C with center A, orthogonal to C.
function orthogonalCircle(c, a) {
const ao = sub(a, c);
const h2 = dot(ao, ao);
const r2 = h2 - c[2];
return [...a, r2];
}

// Construction 2.1. Given points A,B ∈ D, construct the hyperbolic geodesic AB.
// Equivalently, given two points A,B and a circle C∞ with center O, construct the unique
// circle through A,B that is orthogonal to C∞
function geodesic(c, a, b) {
if (Math.abs(a[0]*b[1] - b[0]*a[1]) < 1.0e-14) {
const n = normalize([b[1] - a[1], a[0] - b[0]]);
return [...n, 0];
}
const ai = circleInvert(c, a);
return circleConstruct3Points(a, b, ai);
}

// invert point about geodesic (either arc or line)
function reflectPG(p, c) {
if (c[2] == 0) {
} else {
return circleInvert(c, p);
}
}

function compass(center, p) {
return [...center, dist_sqr(center, p)];
}

function bisector(a, b) {
const n = sub(b, a);
return [...n, -dot(n, midpoint(a, b))];
}

// Construction of center polygon and reflecting adjacent polygons using
// David E. Joyce's code (https://mathcs.clarku.edu/~djoyce/poincare/).

function constructCenterPolygonVertices(n, k) {
// Initialize P as the center polygon in an n-k regular tiling.
// Let ABC be a triangle in a regular (n,k0-tiling, where
//    A is the center of an n-gon (also center of the disk),
//    B is a vertex of the n-gon, and
//    C is the midpoint of a side of the n-gon adjacent to B.
const angleA = Math.PI/n;
const angleB = Math.PI/k;
const angleC = Math.PI/2.0;
// For a regular tiling, we need to compute the distance s from A to B.
const sinA = Math.sin(angleA);
const sinB = Math.sin(angleB);
let   s = Math.sin(angleC - angleB - angleA) / Math.sqrt(1.0 - sinB*sinB - sinA*sinA);
// Now determine the coordinates of the n vertices of the n-gon.
// They're all at distance s from the center of the Poincare disk.
const P = [];
for (let i=0; i<n; ++i) {
P[i] = [s * Math.cos((3+2*i)*angleA),s * Math.sin((3+2*i)*angleA)];
}

if (Math.abs(centerOffset) > 0) {
const offset = [centerOffset * Math.cos(3*angleA) * s, centerOffset * Math.sin(3*angleA) * s];
const d = dot(offset, offset);
const invCtr = scale(offset, 1/d);
for (let i=0; i<n; ++i) {
const t = (1/d - 1)/dist_sqr(P[i], invCtr);
P[i] = add(scale(P[i], t), scale(invCtr, (1 - t)));
P[i][0] = -P[i][0]; // Keep chirality. MLA does this.
}
}

return P;
}

function* PoincareDisk(n, k) {
let totalPolygons = 1;

let a = n*(k-3); // polygons in first layer joined by a vertex
let b = n;       // polygons in first layer joined by an edge
let next_a, next_b;

for (let l=1; l<recursion; ++l) {
if (k == 3) {
next_a = a + b;
next_b = (n-6)*a + (n-5)*b;
} else {
next_a = ((n-2)*(k-3) - 1) * a
+ ((n-3)*(k-3) - 1) * b;
next_b = (n-2)*a + (n-3)*b;
}
totalPolygons += a + b;
a = next_a;
b = next_b;
}

// reflect P thru the point or the side indicated by the side s
//  to produce the resulting polygon
const createNextPolygonVertices = (p, s) => {
const polygon = [];
const g = geodesic([0,0,1], p[s], p[(s+1)%n]);
for (let i=0; i<n; ++i) {
const j = (n+s-i+1) % n;
polygon[j] = reflectPG(p[i], g);
}
return polygon;
}

// rule codes
//   0:  initial polygon.  Needs neighbors on all n sides
//   1:  polygon already has 2 neighbors, but one less around corner needed
//   2:  polygon already has 2 neighbors
//   3:  polygon already has 3 neighbors
//   4:  polygon already has 4 neighbors
const P = [{
vertices: constructCenterPolygonVertices(n, k),
rule: 0,
alternating: 0
}];
yield P[0];

for (let j=1, i=0; i<totalPolygons; ++i) {
let r = P[i].rule;
const special = (r==1);
if (special) {
r = 2;
}
const start = (r==4)? 3 : 2;
const quantity = (k==3 && r!=0) ? n-r-1 : n-r;

for (let s=start; s<start+quantity; ++s) {
// Create a polygon adjacent to P[i]
yield P[j] = {
vertices: createNextPolygonVertices(P[i].vertices, s%n),
rule: (k==3 && s==start && r!=0) ? 4 : 3,
alternating: (P[i].alternating == P[0].alternating) ? 1 : 0
};
j++;

let m = special ? 2 : (s==2 && r!=0) ? 1 : 0;
for ( ; m<k-3; ++m) {
// Create a polygon adjacent to P[j-1]
yield P[j] = {
vertices: createNextPolygonVertices(P[j-1].vertices, 1),
rule: (n==3 && m==k-4)? 1 : 2,
alternating: (P[j-1].alternating == P[0].alternating) ? 1 : 0
};
j++;
}
}
}
}

//
// Draw Poincaré disk model
//

// outline
turtle.jump(0, -size);
turtle.circle(size);

const diskIterator = PoincareDisk(n, k);

function walk(i) {
const d = diskIterator.next();
if (!d.done) {
const p = d.value.vertices;
for (let i=0; i<p.length; i++) {
const v0 = p[i];
const v1 = p[(i+1) % p.length];
if (drawFullGeodesics) {
drawGeodesic(geodesic([0,0,1], v0, v1));
} else {
drawGeodesicSegment(v0, v1, geodesic([0,0,1], v0, v1));
}
}
}
return !d.done;
}

//
// Functions to draw geodesics
//

const geodesicsSegmentsDrawn = {};

function drawGeodesicSegment(p1, p2, c) {
// make sure every segment is drawn once
const hash = toFixed([...p1, ...p2]);
if (geodesicsSegmentsDrawn[hash]) {
return;
}
geodesicsDrawn[hash] = true;

drawGeodesicSegmentRec(p1, p2, c);
}

function drawGeodesicSegmentRec(p1, p2, c, rec=0) {
if (rec==0) {
turtle.jump(scale(p1,size));
}
if (dist(p1, p2) > .5/size && c[2] != 0 && rec < 50) {
const m = add(scale(normalize(sub(midpoint(p1, p2), c)), Math.sqrt(c[2])), c);
drawGeodesicSegmentRec(p1, m, c, rec+1);
drawGeodesicSegmentRec(m, p2, c, rec+1);
} else {
turtle.goto(scale(p2,size));
}
}

const geodesicsDrawn = {};

function drawGeodesic(c) {
// make sure every geodesic is drawn once

if (c[2] == 0) { c[0] * Math.sign(c[1]); c[1] = Math.abs(c[1]); }
const hash = toFixed(c);
if (geodesicsDrawn[hash]) {
return;
}
geodesicsDrawn[hash] = true;

if (c[2] != 0) {
const r = Math.sqrt(c[2]);
turtle.jump(scale([c[0],c[1]-r],size));
turtle.circle(r*size);
} else {
turtle.jump(scale([c[1], -c[0]], 200));
turtle.goto(scale([c[1], -c[0]], -200));
}
}

function toFixed(n) {
let s = '';
n.forEach(i => {
i = Math.abs(i) < 1e-4 ? 0 : i;
s += `\${i.toFixed(3)}-`;
});
return s;
}

//
// 2D Vector math
//

function cross(a, b) { return [a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]]; }
function midpoint(a, b) { return [(a[0]+b[0])*.5, (a[1]+b[1])*.5]; }
function equal(a,b) { return .001>dist_sqr(a,b); }
function scale(a,b) { return [a[0]*b,a[1]*b]; }
function add(a,b) { return [a[0]+b[0],a[1]+b[1]]; }
function sub(a,b) { return [a[0]-b[0],a[1]-b[1]]; }
function dot(a,b) { return a[0]*b[0]+a[1]*b[1]; }
function dist_sqr(a,b) { return (a[0]-b[0])**2+(a[1]-b[1])**2; }
function dist(a,b) { return Math.sqrt(dist_sqr(a,b)); }
function length(a) { return Math.sqrt(dot(a,a)); }
function normalize(a) { return scale(a, 1/length(a)); }
function lerp(a,b,t) { return [a[0]*(1-t)+b[0]*t,a[1]*(1-t)+b[1]*t]; }```