const size = 95;
const hideBelowRadius = .8; //min=.1 max=10 step=.1
const hideAboveRadius = 100; //min=.1 max=100 step=.1
const generations = 2; //min=1 max=30 step=1
const firstChildRatio = .5; //min=0 max=1 step=.01
const secondChildRatio = 1; //min=0 max=1 step=.01
const rotateGeneration = 0; //min=0 max=1 step=.01
const mutation = 0; //min=0 max=1 step=.01 Setting mutation to 1 makes every sibling in generation to have completely random child ratios

const showMaxR = hideAboveRadius;
const minR = hideBelowRadius;

// You can find the Turtle API reference here: https://turtletoy.net/syntax
Canvas.setpenopacity(-.8);

// Global code will be evaluated once.
init();
const turtle = new Turtle();

class CircleGenerator {
    constructor(firstChildRatio = .5, secondChildRatio = 1, mutation = 0) {
        this.firstChildRatio = firstChildRatio;
        this.secondChildRatio = secondChildRatio;
        this.mutation = mutation;
        this.lastGeneration = 0;
    }
    generate(originalCircle, generation = 0) {
        //for the sake of aesthetics I consider sibblings to have the same mutations for their children
        if(generation != this.lastGeneration) {
            const firstChildMutation = 1 + this.mutation * (Math.random() * 2 - 1);
            const secondChildMutation = 1 + this.mutation * (Math.random() * 2 - 1);
            this.firstChildRatio = N.clamp(this.firstChildRatio * firstChildMutation, 0, 1);
            this.secondChildRatio = N.clamp(this.secondChildRatio * secondChildMutation, 0, 1);
            this.lastGeneration = generation;
        }
        if(this.mutation == 1) {
            this.firstChildRatio = Math.random();
            this.secondChildRatio = Math.random();
        }
        const transpose = [...originalCircle.c];
        const circle = new Circle([0,0], originalCircle.b);

        if(circle.r < minR * 2) return [];

        const rA = circle.r * this.firstChildRatio;
        const xA = rA - circle.r;
        
        const rB = Math.abs(
            (circle.r - rA) * this.secondChildRatio // max of rB times ratio
        );
        
        const circleA = new Circle(V.add(circle.c, [xA, 0]), 1/rA);

        const circleD = new Circle(circleA.c, 1/(rA + rB));
        const circleE = new Circle(circle.c, 1/(circle.r - rB));
        const intersections = Intersection.circles(circleD.c, circleD.r, circleE.c, circleE.r);
        
        const circleB = new Circle( this.secondChildRatio < 0? intersections.point_2: intersections.point_1, 1/rB);
        
        circleA.c = V.add( V.trans(V.rot2d(generation * 2 * Math.PI * rotateGeneration), circleA.c), transpose);
        circleB.c = V.add( V.trans(V.rot2d(generation * 2 * Math.PI * rotateGeneration), circleB.c), transpose);
        
        return [circleA, circleB]
    }
}

const c = new Circle([0,0], 1/size);
const cg = new CircleGenerator(firstChildRatio, secondChildRatio, mutation);
const ag = new ApollonianGasket(turtle, c, cg.generate.bind(cg));

ag.circles[0].draw(turtle, showMaxR);

// The walk function will be called until it returns false.
function walk(i) {
    return ag.step();
}

function ApollonianGasket(turtle, circle, generatorFn) {
    class ApollonianGasket {
        constructor(turtle, circle, generatorFn) {
            this.turtle = turtle;
            this.queues = [ ];
            this.generatorFn = generatorFn;
            
            this.addQueue(circle, 0);
        }
        addCircle(circle, triplet) {
            this.drawCircle(circle);
            this.queues[0][0].push(circle);
            triplet.forEach((e, i, a) => 
                this.queues[0][1].push([e, a[(i+1)%a.length], circle])
            );
        }
        drawCircle(circle) {
            circle.draw(this.turtle, showMaxR);
        }
        addQueue(circle, generation = 0) {
            if(generation > generations - 1) return;
            const c = new Circle([...circle.c], -circle.b);
            const newCircles = this.generatorFn(c, generation);
            if(newCircles) {
                newCircles.forEach(cc => this.drawCircle(cc));
                this.queues.push([
                    [c, ...newCircles],
                    [[c, ...newCircles]],
                    generation
                ]);
            }
        }
        hasNext() {
            return this.queues[0] !== undefined && this.queues[0][1] !== undefined && this.queues[0][1].length > 0;
        }
        step() {
            if(this.queues[0] === undefined) {
                return;
            }
            const triplet = this.queues[0][1].shift();
            
            // Generate new circles based on Descartes' theorem
            for (let candidate of complexDescartes(triplet, descartes(...triplet))) {
                if (candidate.r < minR) continue;
                
                if (validate(candidate, triplet)) {
                    this.addCircle(candidate, triplet);
                }
            }
            
            if(this.queues[0][1].length == 0) {
                this.queues[0][0].forEach((c,i) => {
                    if(i === 0) return;
                    this.addQueue(c, this.queues[0][2] + 1);
                });
                
                this.queues.shift();
            }
            
            return this.hasNext();
        }
        get circles() {
            return this.queues[0][0];
        }
    }
    return new ApollonianGasket(turtle, circle, generatorFn);
}
function validate(candidate, triplet) {
    if(ag.circles.some(c => V.approx(candidate.c, c.c) && N.approx(candidate.r, c.r))) return false;
    // Check if all 4 circles are mutually tangential
    return !triplet.some(c => !isTangent(candidate, c));
}

function isTangent(c1, c2) {
    let d = V.len(V.sub(c1.c, c2.c));
    // Tangency check based on distances and radii
    return N.approx(d, c1.r+c2.r) || N.approx(d, Math.abs(c2.r - c1.r))
}

// Complex calculations based on Descartes' theorem for circle generation
// https://en.wikipedia.org/wiki/Descartes%27_theorem
function complexDescartes(circles, k4) {
    const zk = circles.map(c => C.scale(c.c, c.b));
    const sum = zk.reduce((a,c)=>C.add(a,c),[0,0]);
    const product = zk.map((e, i, a) => C.mult(e, a[(i+1)%a.length])).reduce((a,c) => C.add(a,c), [0,0]);
    const root = C.scale(C.sqrt(product), 2);
    
    return k4.flatMap(k => [C.add, C.sub].map(fn => new Circle(C.scale(fn(sum, root), 1/k), k)));
}

// Calculate curvatures (k-values) for new circles using Descartes' theorem
function descartes(...circles) {
    const sum = circles.reduce((a, c) => a + c.b, 0);
    const product = Math.abs(circles.map((e,i,a) => e.b * a[(i+1)%a.length].b).reduce((a,c) => a + c, 0));
    let root = 2 * product**.5;
    return [sum + root, sum - root];
}

function Circle(location, bend) {
    class Circle {
        constructor(location, bend) {
            this.c = location;
            this.r = Math.abs(1/bend);
            this.b = bend;
        }
        draw(turtle, showMaxR) {
            if(this.r > showMaxR) return;
            PT.drawTour(turtle, PT.circle(this.r).map(pt => V.add(pt, this.c)));
        }
    }
    return new Circle(location, bend);
}

function init() {
    ///////////////////////////////////////////////////////
    // Vector functions - Created by Jurgen Westerhof 2024
    // https://turtletoy.net/turtle/d068ad6040
    ///////////////////////////////////////////////////////
    class Vector {
        static add  (a,b) { return a.map((v,i)=>v+b[i]); }
        static sub  (a,b) { return a.map((v,i)=>v-b[i]); }
        static mul  (a,b) { return a.map((v,i)=>v*b[i]); }
        static div  (a,b) { return a.map((v,i)=>v/b[i]); }
        static scale(a,s) { return a.map(v=>v*s); }
    
        static det(m)                { return m.length == 1? m[0][0]: m.length == 2 ? m[0][0]*m[1][1]-m[0][1]*m[1][0]: m[0].reduce((r,e,i) => r+(-1)**(i+2)*e*this.det(m.slice(1).map(c => c.filter((_,j) => i != j))),0); }
        static angle(a)              { return Math.PI - Math.atan2(a[1], -a[0]); } //compatible with turtletoy heading
        static rot2d(angle)          { return [[Math.cos(angle), -Math.sin(angle)], [Math.sin(angle), Math.cos(angle)]]; }
        static rot3d(yaw,pitch,roll) { return [[Math.cos(yaw)*Math.cos(pitch), Math.cos(yaw)*Math.sin(pitch)*Math.sin(roll)-Math.sin(yaw)*Math.cos(roll), Math.cos(yaw)*Math.sin(pitch)*Math.cos(roll)+Math.sin(yaw)*Math.sin(roll)],[Math.sin(yaw)*Math.cos(pitch), Math.sin(yaw)*Math.sin(pitch)*Math.sin(roll)+Math.cos(yaw)*Math.cos(roll), Math.sin(yaw)*Math.sin(pitch)*Math.cos(roll)-Math.cos(yaw)*Math.sin(roll)],[-Math.sin(pitch), Math.cos(pitch)*Math.sin(roll), Math.cos(pitch)*Math.cos(roll)]]; }
        static trans(matrix,a)       { return a.map((v,i) => a.reduce((acc, cur, ci) => acc + cur * matrix[ci][i], 0)); }
        //Mirror vector a in a ray through [0,0] with direction mirror
        static mirror2d(a,mirror)    { return [Math.atan2(...mirror)].map(angle => this.trans(this.rot2d(angle), this.mul([-1,1], this.trans(this.rot2d(-angle), a)))).pop(); }

        static approx(a,b,p) { return this.len(this.sub(a,b)) < (p === undefined? .001: p); }
        static norm  (a)     { return this.scale(a,1/this.len(a)); }
        static len   (a)     { return Math.hypot(...a); }
        static lenSq (a)     { return a.reduce((a,c)=>a+c**2,0); }
        static lerp  (a,b,t) { return a.map((v, i) => v*(1-t) + b[i]*t); }
        static dist  (a,b)   { return Math.hypot(...this.sub(a,b)); }
        
        static dot  (a,b)   { return a.reduce((a,c,i) => a+c*b[i], 0); }
        static cross(...ab) { return ab[0].map((e, i) => ab.map(v => v.filter((ee, ii) => ii != i))).map((m,i) => (i%2==0?-1:1)*this.det(m)); }
    }
    this.V = Vector;
    
    class Intersection2D {
        //a-start, a-direction, b-start, b-direction
        //returns false on no intersection or [[intersection:x,y], scalar a-direction, scalar b-direction
        static info(as, ad, bs, bd) {
            const d = V.sub(bs, as), det = -V.det([bd, ad]);
            if(det === 0) return false;
            const res = [V.det([d, bd]) / det, V.det([d, ad]) / det];
            return [V.add(as, V.scale(ad, res[0])), ...res];
        }
        static ray(a, b, c, d) {
            return this.info(a, b, c, d);
        }
        static segment(a,b,c,d, inclusiveStart = true, inclusiveEnd = true) {
            const i = this.info(a, V.sub(b, a), c, V.sub(d, c));
            return i === false? false: (
                (inclusiveStart? 0<=i[1] && 0<=i[2]: 0<i[1] && 0<i[2])
             && (inclusiveEnd?   i[1]<=1 && i[2]<=1: i[1]<1 && i[2]<1)
            )?i[0]:false;
        }
        static tour(tour, segmentStart, segmentDirection) {
            return tour.map((e, i, a) => [i, this.info(e, V.sub(a[(i+1)%a.length], e), segmentStart, segmentDirection)])
                .filter(e => e[1] !== false && 0 <= e[1][1] && e[1][1] <= 1)
                .filter(e => 0 <= e[1][2])// && e[1][2] <= 1)
                .map(e => ({
                    position: e[1][0],
                    tourIndex: e[0],
                    tourSegmentPortion: e[1][1],
                    segmentPortion: e[1][2],
                })
            );
        }
        static inside(tour, pt) {
            return tour.map((e,i,a) => this.segment(e, a[(i+1)%a.length], pt, [Number.MAX_SAFE_INTEGER, 0], true, false)).filter(e => e !== false).length % 2 == 1
        }
        static circles(centerA, radiusA, centerB, radiusB) {
            // Start constructing the response object.
            const result = {
                intersect_count: 0,
                intersect_occurs: true,
                one_is_in_other: false,
                are_equal: false,
                point_1: [null, null],
                point_2: [null, null],
            };
        
            // Get vertical and horizontal distances between circles.
            const dx = centerB[0] - centerA[0];
            const dy = centerB[1] - centerA[1];
        
            // Calculate the distance between the circle centers as a straight line.
            const dist = Math.hypot(dy, dx);
        
            // Check if circles intersect.
            if (dist > radiusA + radiusB) {
                result.intersect_occurs = false;
            }
        
            // Check one circle isn't inside the other.
            if (dist < Math.abs(radiusA - radiusB) && !N.approx(dist, Math.abs(radiusA - radiusB))) {
                result.intersect_occurs = false;
                result.one_is_in_other = true;
            }
        
            // Check if circles are the same.
            if (V.approx(centerA, centerB) && radiusA === radiusB) {
                result.are_equal = true;
            }
        
            // Find the intersection points
            if (result.intersect_occurs) {
                // Centroid is the pt where two lines cross. A line between the circle centers
                // and a line between the intersection points.
                const centroid = (radiusA**2 - radiusB**2 + dist * dist) / (2.0 * dist);

                // Get the coordinates of centroid.
                const x2 = centerA[0] + (dx * centroid) / dist;
                const y2 = centerA[1] + (dy * centroid) / dist;
                
                const prec = 10000;
                // Get the distance from centroid to the intersection points.
                const h = (Math.round(radiusA**2 * prec)/prec - Math.round(centroid**2 * prec)/prec)**.5;
        
                // Get the x and y dist of the intersection points from centroid.
                const rx = -dy * (h / dist);
                const ry = dx * (h / dist);
                
                // Get the intersection points.
                result.point_1 = [x2 + rx, y2 + ry];
                result.point_2 = [x2 - rx, y2 - ry];

                // Add intersection count to results
                if (result.are_equal) {
                    result.intersect_count = null;
                } else if (result.point_1.x === result.point_2.x && result.point_1.y === result.point_2.y) {
                    result.intersect_count = 1;
                } else {
                    result.intersect_count = 2;
                }
            }
            return result;
        }
    }
    this.Intersection = Intersection2D;
    
    class PathTools {
        static bezier(p1, cp1, cp2, p2, steps = null) {steps = (steps === null? (V.len(V.sub(cp1, p1)) + V.len(V.sub(cp2, cp1)) + V.len(V.sub(p2, cp2))) | 0: steps) - 1;return Array.from({length: steps + 1}).map((v, i, a, f = i/steps) => [[V.lerp(p1, cp1, f),V.lerp(cp1, cp2, f),V.lerp(cp2, p2, f)]].map(v => V.lerp(V.lerp(v[0], v[1], f), V.lerp(v[1], v[2], f), f))[0]);}
        // https://stackoverflow.com/questions/18655135/divide-bezier-curve-into-two-equal-halves#18681336
        static splitBezier(p1, cp1, cp2, p2, t=.5) {const e = V.lerp(p1, cp1, t);const f = V.lerp(cp1, cp2, t);const g = V.lerp(cp2, p2, t);const h = V.lerp(e, f, t);const j = V.lerp(f, g, t);const k = V.lerp(h, j, t);return [[p1, e, h, k], [k, j, g, p2]];}
        static circular(radius,verticeCount,rotation=0) {return Array.from({length: verticeCount}).map((e,i,a,f=i*2*Math.PI/verticeCount+rotation) => [radius*Math.cos(f),radius*Math.sin(f)])}
        static circle(r){return this.circular(r,Math.max(12, r*2*Math.PI|0));}
        static draw(turtle, path) {path.forEach((pt, i) => turtle[i==0?'jump':'goto'](pt));}
        static drawTour(turtle, path) {this.draw(turtle, path.concat([path[0]]));}
        static drawPoint(turtle, pt, r = .1) {this.drawTour(turtle, this.circle(r).map(e => V.add(e, pt)));}
        static drawArrow(turtle, s, d, width = 6, length = 3) {
            turtle.jump(s);
            const arrowHeadBase = V.add(s,d);
            turtle.goto(arrowHeadBase);
            turtle.goto(V.add(arrowHeadBase, V.trans(V.rot2d(-V.angle(d)), [-length, width/2])));
            turtle.jump(V.add(arrowHeadBase, V.trans(V.rot2d(-V.angle(d)), [-length, -width/2])));
            turtle.goto(arrowHeadBase);
        }
    }
    
    this.PT = PathTools;
    
    class Complex {
        static add(a,b)     { return V.add(a,b); }
        static sub(a,b)     { return V.sub(a,b); }
        static scale(a,s)   { return V.scale(a,s); }
        static mult(a,b)    { return [a[0]*b[0]-a[1]*b[1],a[0]*b[1]+a[1]*b[0]]; }
        static sqrt(a)      { return [[Math.hypot(...a)**.5, Math.atan2(...a.reverse()) / 2]].map(ra => [ra[0]*Math.cos(ra[1]), ra[0]*Math.sin(ra[1])]).pop(); }
    }
    this.C = Complex;
    
    class Numbers {
        static approx(a,b,p)        { return Math.abs(a-b) < (p === undefined? .001: p); }
        static clamp(a, min, max)   { return Math.min(Math.max(a, min), max); }
    }
    this.N = Numbers;
}