// Aperiodic Hat Monotiles. Created by Reinder Nijhoff 2023
// @reindernijhoff
//
// https://turtletoy.net/turtle/61a52c764a
// 
// Code adapted from: https://github.com/isohedral/hatviz
// BSD 3-Clause License
// Copyright (c) 2023, Craig S. Kaplan
//

const turtle = new Turtle();

const levels = 3; // min=1, max=6, step=1
const scale = 6; // min=0.5, max=10, step=0.001
const draw  = 1; // min=1, max=3, step=1 (Draw Hats, Draw Supertiles, Draw All)
const hatching = 1; // min=0, max=1, step=1 (No, Yes)

const polygons = new Polygons();
const iterator = createEinsteinTiling(levels, draw, scale);

function walk(i) {
	const d = iterator.next();
	if (!d.done) {
	    const v = d.value;
	    
        const poly = polygons.create();
        poly.addPoints(...v.polygon);
        poly.addOutline();
        if (hatching && v.geom.children.length === 0) {
            switch(v.geom.fill) {
                case 0: poly.addHatching( .8,1.0);
                        break;
                case 1: poly.addHatching(-.8,0.7);
                        break;
                case 2: poly.addHatching( .1,0.3);
                        break;
            }
        }
        polygons.draw(turtle, poly, false);
        
		return true;
	} else {
		return false;
	}
}

//
// An aperiodic monotile
// David Smith, Joseph Samuel Myers, Craig S. Kaplan, and Chaim Goodman-Strauss, 2023
// https://cs.uwaterloo.ca/~csk/hat/
// 
// Code adapted from: https://github.com/isohedral/hatviz
// BSD 3-Clause License
// Copyright (c) 2023, Craig S. Kaplan
//

function* createEinsteinTiling(levels, draw = 1, scale = 1) {
	const r3 = Math.sqrt(3);
	const hr3 = r3 / 2;
	const ident = [1, 0, 0, 0, 1, 0];

    const scl2   = (a,b)   => [a[0]*b, a[1]*b];
    const add2   = (a,b)   => [a[0]+b[0], a[1]+b[1]];
    const sub2   = (a,b)   => [a[0]-b[0], a[1]-b[1]];
    const rot2   = (a)     => [Math.cos(a), -Math.sin(a), Math.sin(a), Math.cos(a)];
	const hexPt  = (x, y)  => [x + 0.5 * y, hr3 * y];
	const intersect = (p1, q1, p2, q2) => {
		const d = (q2[1] - p2[1]) * (q1[0] - p1[0]) - (q2[0] - p2[0]) * (q1[1] - p1[1]);
		const uA = ((q2[0] - p2[0]) * (p1[1] - p2[1]) - (q2[1] - p2[1]) * (p1[0] - p2[0])) / d;
		return [p1[0] + uA * (q1[0] - p1[0]), p1[1] + uA * (q1[1] - p1[1])];
	}
    
	// Affine matrix inverse
    const inv = (T) => {
		const det = T[0] * T[4] - T[1] * T[3];
		return [ T[4] / det, -T[1] / det, (T[1] * T[5] - T[2] * T[4]) / det,
			     -T[3] / det, T[0] / det, (T[2] * T[3] - T[0] * T[5]) / det ];
	}
	// Affine matrix multiply
	const mul = (A, B) => [
		    A[0] * B[0] + A[1] * B[3],
			A[0] * B[1] + A[1] * B[4],
			A[0] * B[2] + A[1] * B[5] + A[2],

			A[3] * B[0] + A[4] * B[3],
			A[3] * B[1] + A[4] * B[4],
			A[3] * B[2] + A[4] * B[5] + A[5]
		];
	// Rotation matrix
	const trot = (a) => [Math.cos(a), -Math.sin(a), 0, Math.sin(a), Math.cos(a), 0];
	// Translation matrix
	const ttrans = (p) => [1, 0, p[0], 0, 1, p[1]];
    const rotAbout = (p, ang) => mul(ttrans(p), mul(trot(ang), ttrans([-p[0], -p[1]])));
	// Matrix * point
	const transPt = (M, P) => [M[0] * P[0] + M[1] * P[1] + M[2], M[3] * P[0] + M[4] * P[1] + M[5]];
	// Match unit interval to line segment p->q
	const matchSeg = (p, q) => [q[0] - p[0], p[1] - q[1], p[0], q[1] - p[1], q[0] - p[0], p[1]];
	// Match line segment p1->q1 to line segment p2->q2
	const matchTwo = (p1, q1, p2, q2) => mul(matchSeg(p2, q2), inv(matchSeg(p1, q1)));


	class Geom {
		constructor(pgon, fill = 0) {
			this.shape = pgon;
			this.fill = fill
			this.width = 1.0;
			this.children = [];
		}

		addChild(T, geom) {
			this.children.push({
				T: T,
				geom: geom
			});
		}

		evalChild(n, i) {
			return transPt(this.children[n].T, this.children[n].geom.shape[i]);
		}

		recentre() {
			let tr = [0, 0];
			this.shape.forEach( p => tr = add2(tr, p));
			tr = scl2( tr, -1 / this.shape.length);

			this.shape = this.shape.map (p => add2(p, tr));

			const M = ttrans(tr);
			for (let ch of this.children) {
				ch.T = mul(M, ch.T);
			}
		}
	}
	

	const hat_outline = [
		hexPt(0, 0), hexPt(-1, -1), hexPt(0, -2), hexPt(2, -2),
		hexPt(2, -1), hexPt(4, -2), hexPt(5, -1), hexPt(4, 0),
		hexPt(3, 0), hexPt(2, 2), hexPt(0, 3), hexPt(0, 2),
		hexPt(-1, 2)
	];

	const H1_hat = new Geom(hat_outline, 0);
	const H_hat = new Geom(hat_outline, 1);
	const T_hat = new Geom(hat_outline, 2);
	const P_hat = new Geom(hat_outline, 3);
	const F_hat = new Geom(hat_outline, 4);

	function constructPatch(H, T, P, F) {
		const rules = [
			['H'],
			[0, 0, 'P', 2],
			[1, 0, 'H', 2],
			[2, 0, 'P', 2],
			[3, 0, 'H', 2],
			[4, 4, 'P', 2],
			[0, 4, 'F', 3],
			[2, 4, 'F', 3],
			[4, 1, 3, 2, 'F', 0],
			[8, 3, 'H', 0],
			[9, 2, 'P', 0],
			[10, 2, 'H', 0],
			[11, 4, 'P', 2],
			[12, 0, 'H', 2],
			[13, 0, 'F', 3],
			[14, 2, 'F', 1],
			[15, 3, 'H', 4],
			[8, 2, 'F', 1],
			[17, 3, 'H', 0],
			[18, 2, 'P', 0],
			[19, 2, 'H', 2],
			[20, 4, 'F', 3],
			[20, 0, 'P', 2],
			[22, 0, 'H', 2],
			[23, 4, 'F', 3],
			[23, 0, 'F', 3],
			[16, 0, 'P', 2],
			[9, 4, 0, 2, 'T', 2],
			[4, 0, 'F', 3]
		];

		ret = new Geom([], null, null);
		ret.width = H.width;
		shapes = {
			'H': H,
			'T': T,
			'P': P,
			'F': F
		};

		for (let r of rules) {
			if (r.length == 1) {
				ret.addChild(ident, shapes[r[0]]);
			} else if (r.length == 4) {
				const poly = ret.children[r[0]].geom.shape;
				const T = ret.children[r[0]].T;
				const P = transPt(T, poly[(r[1] + 1) % poly.length]);
				const Q = transPt(T, poly[r[1]]);
				const nshp = shapes[r[2]];
				const npoly = nshp.shape;

				ret.addChild(matchTwo(npoly[r[3]], npoly[(r[3] + 1) % npoly.length], P, Q), nshp);
			} else {
				const chP = ret.children[r[0]];
				const chQ = ret.children[r[2]];

				const P = transPt(chQ.T, chQ.geom.shape[r[3]]);
				const Q = transPt(chP.T, chP.geom.shape[r[1]]);
				const nshp = shapes[r[4]];
				const npoly = nshp.shape;

				ret.addChild(matchTwo(npoly[r[5]], npoly[(r[5] + 1) % npoly.length], P, Q), nshp);
			}
		}

		return ret;
	}

	function constructMetatiles(patch) {
		const bps1 = patch.evalChild(8, 2);
		const bps2 = patch.evalChild(21, 2);
		const rbps = transPt(rotAbout(bps1, -2.0 * Math.PI / 3.0), bps2);

		const p72 = patch.evalChild(7, 2);
		const p252 = patch.evalChild(25, 2);

		const llc = intersect(bps1, rbps, patch.evalChild(6, 2), p72);
		let w = sub2(patch.evalChild(6, 2), llc);

		const new_H_outline = [llc, bps1];
		w = transPt(trot(-Math.PI / 3), w);
		new_H_outline.push(add2(new_H_outline[1], w));
		new_H_outline.push(patch.evalChild(14, 2));
		w = transPt(trot(-Math.PI / 3), w);
		new_H_outline.push(sub2(new_H_outline[3], w));
		new_H_outline.push(patch.evalChild(6, 2));

		const new_H = new Geom(new_H_outline);
		new_H.width = patch.width * 2;
		for (let ch of [0, 9, 16, 27, 26, 6, 1, 8, 10, 15]) {
			new_H.addChild(patch.children[ch].T, patch.children[ch].geom);
		}

		const new_P_outline = [p72, add2(p72, sub2(bps1, llc)), bps1, llc];
		const new_P = new Geom(new_P_outline);
		new_P.width = patch.width * 2;
		for (let ch of [7, 2, 3, 4, 28]) {
			new_P.addChild(patch.children[ch].T, patch.children[ch].geom);
		}

		const new_F_outline = [bps2, patch.evalChild(24, 2), patch.evalChild(25, 0), p252, add2(p252, sub2(llc, bps1))];
		const new_F = new Geom(new_F_outline);
		new_F.width = patch.width * 2;
		for (let ch of [21, 20, 22, 23, 24, 25]) {
			new_F.addChild(patch.children[ch].T, patch.children[ch].geom);
		}

		const AAA = new_H_outline[2];
		const BBB = add2(new_H_outline[1], sub2(new_H_outline[4], new_H_outline[5]));
		const CCC = transPt(rotAbout(BBB, -Math.PI / 3), AAA);
		const new_T_outline = [BBB, CCC, AAA];
		const new_T = new Geom(new_T_outline);
		new_T.width = patch.width * 2;
		new_T.addChild(patch.children[11].T, patch.children[11].geom);

		new_H.recentre();
		new_P.recentre();
		new_F.recentre();
		new_T.recentre();

		return [new_H, new_T, new_P, new_F]
	}

	// init start tiles
	const H_outline = [[0, 0], [4, 0], [4.5, hr3], [2.5, 5 * hr3], [1.5, 5 * hr3], [-0.5, hr3]];
	const H_init = new Geom(H_outline);
	H_init.width = 2;
	H_init.addChild(matchTwo(hat_outline[5], hat_outline[7], H_outline[5], H_outline[0]), H_hat);
	H_init.addChild(matchTwo(hat_outline[9], hat_outline[11], H_outline[1], H_outline[2]), H_hat);
	H_init.addChild(matchTwo(hat_outline[5], hat_outline[7], H_outline[3], H_outline[4]), H_hat);
	H_init.addChild(mul(ttrans([2.5, hr3]), mul([-0.5, -hr3, 0, hr3, -0.5, 0], [0.5, 0, 0, 0, -0.5, 0])), H1_hat);

	const T_outline = [[0, 0], [3, 0], [1.5, 3 * hr3]];
	const T_init = new Geom(T_outline);
	T_init.width = 2;
	T_init.addChild([0.5, 0, 0.5, 0, 0.5, hr3], T_hat);

	const P_outline = [[0, 0], [4, 0], [3, 2 * hr3], [-1, 2 * hr3]];
	const P_init = new Geom(P_outline);
	P_init.width = 2;
	P_init.addChild([0.5, 0, 1.5, 0, 0.5, hr3], P_hat);
	P_init.addChild(mul(ttrans([0, 2 * hr3]), mul([0.5, hr3, 0, -hr3, 0.5, 0], [0.5, 0.0, 0.0, 0.0, 0.5, 0.0])), P_hat);

	const F_outline = [[0, 0], [3, 0], [3.5, hr3], [3, 2 * hr3], [-1, 2 * hr3]];
	const F_init = new Geom(F_outline);
	F_init.width = 2;
	F_init.addChild([0.5, 0, 1.5, 0, 0.5, hr3], F_hat);
	F_init.addChild(mul(ttrans([0, 2 * hr3]), mul([0.5, hr3, 0, -hr3, 0.5, 0], [0.5, 0.0, 0.0, 0.0, 0.5, 0.0])), F_hat);

	let tiles = [H_init, T_init, P_init, F_init];

	// create all tiles for level
	for (let i = 0; i < levels; i++) {
		const patch = constructPatch(...tiles);
		tiles = constructMetatiles(patch);
	}

	// generator code
	const queue = [{
		T: [scale, 0, 0, 0, scale, 0],
		geom: tiles[0],
		level: levels
	}];
	do {
		const t = queue.pop();
		if (t.level >= 0) {
			for (let g of t.geom.children) {
				queue.push({
					T: mul(t.T, g.T),
					geom: g.geom,
					level: t.level - 1
				});
			}
		} 
		if ( (t.level < 0 && draw == 1) || (t.level >= 0 && draw == 2) || draw == 3) {
    		yield {
    		    geom: t.geom,
    		    polygon: t.geom.shape.map(p => transPt(t.T, p))
    		};
		}
	} while (queue.length > 0);
}

////////////////////////////////////////////////////////////////
// Polygon Clipping utility code - Created by Reinder Nijhoff 2019
// (Polygon binning by Lionel Lemarie 2021)
// https://turtletoy.net/turtle/a5befa1f8d
////////////////////////////////////////////////////////////////
function Polygons(){const t=[],s=25,e=Array.from({length:s**2},t=>[]),n=class{constructor(){this.cp=[],this.dp=[],this.aabb=[]}addPoints(...t){let s=1e5,e=-1e5,n=1e5,h=-1e5;(this.cp=[...this.cp,...t]).forEach(t=>{s=Math.min(s,t[0]),e=Math.max(e,t[0]),n=Math.min(n,t[1]),h=Math.max(h,t[1])}),this.aabb=[s,n,e,h]}addSegments(...t){t.forEach(t=>this.dp.push(t))}addOutline(){for(let t=0,s=this.cp.length;t<s;t++)this.dp.push(this.cp[t],this.cp[(t+1)%s])}draw(t){for(let s=0,e=this.dp.length;s<e;s+=2)t.jump(this.dp[s]),t.goto(this.dp[s+1])}addHatching(t,s){const e=new n;e.cp.push([-1e5,-1e5],[1e5,-1e5],[1e5,1e5],[-1e5,1e5]);const h=Math.sin(t)*s,o=Math.cos(t)*s,a=200*Math.sin(t),i=200*Math.cos(t);for(let t=.5;t<150/s;t++)e.dp.push([h*t+i,o*t-a],[h*t-i,o*t+a]),e.dp.push([-h*t+i,-o*t-a],[-h*t-i,-o*t+a]);e.boolean(this,!1),this.dp=[...this.dp,...e.dp]}inside(t){let s=0;for(let e=0,n=this.cp.length;e<n;e++)this.segment_intersect(t,[.1,-1e3],this.cp[e],this.cp[(e+1)%n])&&s++;return 1&s}boolean(t,s=!0){const e=[];for(let n=0,h=this.dp.length;n<h;n+=2){const h=this.dp[n],o=this.dp[n+1],a=[];for(let s=0,e=t.cp.length;s<e;s++){const n=this.segment_intersect(h,o,t.cp[s],t.cp[(s+1)%e]);!1!==n&&a.push(n)}if(0===a.length)s===!t.inside(h)&&e.push(h,o);else{a.push(h,o);const n=o[0]-h[0],i=o[1]-h[1];a.sort((t,s)=>(t[0]-h[0])*n+(t[1]-h[1])*i-(s[0]-h[0])*n-(s[1]-h[1])*i);for(let n=0;n<a.length-1;n++)(a[n][0]-a[n+1][0])**2+(a[n][1]-a[n+1][1])**2>=.001&&s===!t.inside([(a[n][0]+a[n+1][0])/2,(a[n][1]+a[n+1][1])/2])&&e.push(a[n],a[n+1])}}return(this.dp=e).length>0}segment_intersect(t,s,e,n){const h=(n[1]-e[1])*(s[0]-t[0])-(n[0]-e[0])*(s[1]-t[1]);if(0===h)return!1;const o=((n[0]-e[0])*(t[1]-e[1])-(n[1]-e[1])*(t[0]-e[0]))/h,a=((s[0]-t[0])*(t[1]-e[1])-(s[1]-t[1])*(t[0]-e[0]))/h;return o>=0&&o<=1&&a>=0&&a<=1&&[t[0]+o*(s[0]-t[0]),t[1]+o*(s[1]-t[1])]}};return{list:()=>t,create:()=>new n,draw:(n,h,o=!0)=>{reducedPolygonList=function(n){const h={},o=200/s;for(var a=0;a<s;a++){const c=a*o-100,r=[0,c,200,c+o];if(!(n[3]<r[1]||n[1]>r[3]))for(var i=0;i<s;i++){const c=i*o-100;r[0]=c,r[2]=c+o,n[0]>r[2]||n[2]<r[0]||e[i+a*s].forEach(s=>{const e=t[s];n[3]<e.aabb[1]||n[1]>e.aabb[3]||n[0]>e.aabb[2]||n[2]<e.aabb[0]||(h[s]=1)})}}return Array.from(Object.keys(h),s=>t[s])}(h.aabb);for(let t=0;t<reducedPolygonList.length&&h.boolean(reducedPolygonList[t]);t++);h.draw(n),o&&function(n){t.push(n);const h=t.length-1,o=200/s;e.forEach((t,e)=>{const a=e%s*o-100,i=(e/s|0)*o-100,c=[a,i,a+o,i+o];c[3]<n.aabb[1]||c[1]>n.aabb[3]||c[0]>n.aabb[2]||c[2]<n.aabb[0]||t.push(h)})}(h)}}}