// Poincaré disk model #2. Created by Reinder Nijhoff 2021 - @reindernijhoff
//
// https://turtletoy.net/turtle/2341de990a
//
// 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/
// [3] https://www.shadertoy.com/view/3tBGDD
//

Canvas.setpenopacity(.5);

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

if ( n==3 )      { k = Math.max(k,8); }
else if (n == 4) { k = Math.max(k,6); }
else             { k = Math.max(k,4); }

// 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);
}

// Construction 3.1. Given A,B ∈ D, construct the perpendicular bisector of the segment AB. 
// Similarly, construct the midpoint of the segment AB.
function hyperbolicBisector(q, p) {
    const pp = dot(p, p);
    const qq = dot(q, q);
    
    if (pp < 1e-14) { 
        const h2 = 1/qq;
        return [...scale(q,h2), (h2 - 1)];
    } else if (qq < 1e-14) { 
        const h2 = 1/pp;
        return [...scale(p,h2), (h2 - 1.)];
    } else if (Math.abs(pp-qq) < 1e-14) {
        return [...normalize(sub(p,q)), 0];
    }
    
    const d = sub(q, p);
    const k = (1 - dot(p,p))/(dot(d,d) + 2*dot(p,d));
    
    return orthogonalCircle([0,0,1], add(p, scale(d, k)));
}

// invert point about geodesic (either arc or line)
function reflectPG(p, c) {
    if (c[2] == 0) {
        return add(p, scale(c, -2*dot(p, c)));
    } 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;
        const ps = new Polygons();
        const po = ps.create();
        for (let i=0; i<p.length; i++) {
            const v0 = p[i];
            const v1 = p[(i+1) % p.length];
            addGeodesicSegment(po, v0, v1, geodesic([0,0,1], v0, v1));
            if (drawBisectors) {
                drawGeodesic(hyperbolicBisector(v0, v1));
            }
        }
        po.addOutline();
        if (k % 2 == 0) {
            if (d.value.alternating) {
                po.addHatching(Math.PI/4, .5);
                ps.draw(turtle, po);
            }
        } else {
            po.addHatching(Math.PI/4, .1 + Math.random());
            ps.draw(turtle, po);
        }
    }
    return !d.done;
}

function addGeodesicSegment(p, p1, p2, c, rec=0) {
    if (rec==0) {
        p.addPoints(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);
        addGeodesicSegment(p, p1, m, c, rec+1);
        addGeodesicSegment(p, m, p2, c, rec+1);
    } else {
        p.addPoints(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]; }


////////////////////////////////////////////////////////////////
// Polygon Clipping utility code - Created by Reinder Nijhoff 2019
// https://turtletoy.net/turtle/a5befa1f8d
////////////////////////////////////////////////////////////////

function Polygons() {
	const polygonList = [];
	const Polygon = class {
		constructor() {
			this.cp = [];       // clip path: array of [x,y] pairs
			this.dp = [];       // 2d lines [x0,y0],[x1,y1] to draw
			this.aabb = [];     // AABB bounding box
		}
		addPoints(...points) {
		    // add point to clip path and update bounding box
		    let xmin = 1e5, xmax = -1e5, ymin = 1e5, ymax = -1e5;
			(this.cp = [...this.cp, ...points]).forEach( p => {
				xmin = Math.min(xmin, p[0]), xmax = Math.max(xmax, p[0]);
				ymin = Math.min(ymin, p[1]), ymax = Math.max(ymax, p[1]);
			});
		    this.aabb = [(xmin+xmax)/2, (ymin+ymax)/2, (xmax-xmin)/2, (ymax-ymin)/2];
		}
		addSegments(...points) {
		    // add segments (each a pair of points)
		    points.forEach(p => this.dp.push(p));
		}
		addOutline() {
			for (let i = 0, l = this.cp.length; i < l; i++) {
				this.dp.push(this.cp[i], this.cp[(i + 1) % l]);
			}
		}
		draw(t) {
			for (let i = 0, l = this.dp.length; i < l; i+=2) {
				t.jump(this.dp[i]), t.goto(this.dp[i + 1]);
			}
		}
		addHatching(a, d) {
			const tp = new Polygon();
			tp.cp.push([-1e5,-1e5],[1e5,-1e5],[1e5,1e5],[-1e5,1e5]);
			const dx = Math.sin(a) * d,   dy = Math.cos(a) * d;
			const cx = Math.sin(a) * 200, cy = Math.cos(a) * 200;
			for (let i = 0.5; i < 150 / d; i++) {
				tp.dp.push([dx * i + cy,   dy * i - cx], [dx * i - cy,   dy * i + cx]);
				tp.dp.push([-dx * i + cy, -dy * i - cx], [-dx * i - cy, -dy * i + cx]);
			}
			tp.boolean(this, false);
			this.dp = [...this.dp, ...tp.dp];
		}
		inside(p) {
			let int = 0; // find number of i ntersection points from p to far away
			for (let i = 0, l = this.cp.length; i < l; i++) {
				if (this.segment_intersect(p, [0.1, -1000], this.cp[i], this.cp[(i + 1) % l])) {
					int++;
				}
			}
			return int & 1; // if even your outside
		}
		boolean(p, diff = true) {
		    // bouding box optimization by ge1doot.
		    if (Math.abs(this.aabb[0] - p.aabb[0]) - (p.aabb[2] + this.aabb[2]) >= 0 &&
				Math.abs(this.aabb[1] - p.aabb[1]) - (p.aabb[3] + this.aabb[3]) >= 0) return this.dp.length > 0;
				
			// polygon diff algorithm (narrow phase)
			const ndp = [];
			for (let i = 0, l = this.dp.length; i < l; i+=2) {
				const ls0 = this.dp[i];
				const ls1 = this.dp[i + 1];
				// find all intersections with clip path
				const int = [];
				for (let j = 0, cl = p.cp.length; j < cl; j++) {
					const pint = this.segment_intersect(ls0, ls1, p.cp[j], p.cp[(j + 1) % cl]);
					if (pint !== false) {
						int.push(pint);
					}
				}
				if (int.length === 0) {
					// 0 intersections, inside or outside?
					if (diff === !p.inside(ls0)) {
						ndp.push(ls0, ls1);
					}
				} else {
					int.push(ls0, ls1);
					// order intersection points on line ls.p1 to ls.p2
					const cmpx = ls1[0] - ls0[0];
					const cmpy = ls1[1] - ls0[1];
					int.sort( (a,b) =>  (a[0] - ls0[0]) * cmpx + (a[1] - ls0[1]) * cmpy - 
					                    (b[0] - ls0[0]) * cmpx - (b[1] - ls0[1]) * cmpy);
					 
					for (let j = 0; j < int.length - 1; j++) {
						if ((int[j][0] - int[j+1][0])**2 + (int[j][1] - int[j+1][1])**2 >= 0.001) {
							if (diff === !p.inside([(int[j][0]+int[j+1][0])/2,(int[j][1]+int[j+1][1])/2])) {
								ndp.push(int[j], int[j+1]);
							}
						}
					}
				}
			}
			return (this.dp = ndp).length > 0;
		}
		//port of http://paulbourke.net/geometry/pointlineplane/Helpers.cs
		segment_intersect(l1p1, l1p2, l2p1, l2p2) {
			const d   = (l2p2[1] - l2p1[1]) * (l1p2[0] - l1p1[0]) - (l2p2[0] - l2p1[0]) * (l1p2[1] - l1p1[1]);
			if (d === 0) return false;
			const n_a = (l2p2[0] - l2p1[0]) * (l1p1[1] - l2p1[1]) - (l2p2[1] - l2p1[1]) * (l1p1[0] - l2p1[0]);
			const n_b = (l1p2[0] - l1p1[0]) * (l1p1[1] - l2p1[1]) - (l1p2[1] - l1p1[1]) * (l1p1[0] - l2p1[0]);
			const ua = n_a / d;
			const ub = n_b / d;
			if (ua >= 0 && ua <= 1 && ub >= 0 && ub <= 1) {
				return [l1p1[0] + ua * (l1p2[0] - l1p1[0]), l1p1[1] + ua * (l1p2[1] - l1p1[1])];
			}
			return false;
		}
	};
	return {
		list: () => polygonList,
		create: () => new Polygon(),
		draw: (turtle, p, addToVisList=true) => {
			for (let j = 0; j < polygonList.length && p.boolean(polygonList[j]); j++);
			p.draw(turtle);
			if (addToVisList) polygonList.push(p);
		}
	};
}