// equal radii is fairly easy,
// but to handle different sizes,
// must find where 2 spheres 
// intersect each other.
// Trickier!

// see my comment at http://shadertoy.com/view/DdsXW4
void bubbleIntersector(float d, float r0, float rd, out float p, out float r)
{
    // handle equal radii planar case, treat as kinda big sphere (but not too big)
    if (abs(r0 - rd) < 1e-9) { r = 1e3*r0; p = .5*d+r; return; }
    p = d * r0 / (r0 - rd);
    r = sqrt((p - d + (rd*rd - r0*r0) / d)*p + r0*r0);
}
// saved the original r0 > rd version
// TODO handle r0 <= rd
/*
    float rm = max(r0, rd); rd = min(r0, rd); r0 = rm;
    p = d * r0 / (r0 - rd);
    r = sqrt((((p) - d + -abs(rd*rd - r0*r0) / d)*(p) + rm*rm));
*/
// I have most of the nastiness isolated to this one function now.
void bubbleIntersectorRobust(float bd, vec2 p0, vec2 p1, float r0, float r1, out vec2 p2, out float r2)
{
    //float bd = distance(p0, p1);
    // handle equal radii planar case, treat as kinda big sphere (but not too big)
    if (abs(r0 - r1) < 1e-9) { r2 = 1e3*r0; p2 = 1e3 * (p1 - p0) + p0; return; } //d2 = .5*bd+r2; 
    float s = 1.; 
    if (r1 > r0) { s = -1.; float tr = r0; r0 = r1; r1 = tr; vec2 tp = p0; p0 = p1; p1 = tp; } // HACK ensures r0 >= r1
    float d2 = bd * r0 / (r0 - r1);
    r2 = sqrt((d2 - bd + (r1*r1 - r0*r0) / bd)*d2 + r0*r0);
    //r2 = abs(r2); // but it's always positive anyway
    // but the outer loop in mainImage already verified they touch, so the intersection
    // will exist, so r2 will exist (but might be infinite)
    //r2 = min(r2, 1e7); // finite so we don't have to also handle linear case - but already did handle that!
    r2 *= s;
    p2 = d2 / bd * (p1 - p0) + p0;
}
// so now I split off a new one for experimentation
void bubbleIntersectorAttempt(float bd, vec2 p0, vec2 p1, float r0, float r1, out vec2 p2, out float r2)
{
    if (abs(r0 - r1) < 1e-9) { r2 = 1e3*r0; p2 = 1e3 * (p1 - p0) + p0; return; }
    float s = 1.;
    if (r1 > r0) { s = -1.; float tr = r0; r0 = r1; r1 = tr; vec2 tp = p0; p0 = p1; p1 = tp; } // HACK ensures r0 >= r1
    float d2 = bd * r0 / (r0 - r1);
    r2 = sqrt((d2 - bd + (r1*r1 - r0*r0) / bd)*d2 + r0*r0) * s;
    p2 = d2 / bd * (p1 - p0) + p0;
}
// so those output position and signed radius of dividing bubble, handling the inside-out sphere for negative radius is important
// but FIXME I still dislike having to pull the sign back out, but damn the interface is complicated enough already!
// FIXME it gets really confusing trying to keep track of the various ways I refer to the bubbles.  By array index, by P0,P1,P2, R0,Rd,r ...
// and is different for various parts of the code.  I wrote it, and I get really confused by it.
// So once I straighten that out, I hope to go back and redo the identifiers and such.

// simply return signed distance to the divider
float dbubblesplit(vec2 q, float bd, vec2 p0, vec2 p1, float r0, float r1)
{
    float r2; vec2 p2;
    bubbleIntersectorAttempt/*Robust*/(bd, p0, p1, r0, r1, p2, r2);
    float s = sign(r2);
    return (distance(q, p2) - r2 * s) * s; // just return third bubble signed distance
} // so negative result means not actually inside this bubble b0, but in other b1

#define IZERO min(0, iFrame)

// Brute force is not the way.
// In a real situation, would likely set
// N to be the maximum neighbor count.
// That at least prevents N from being large,
// which makes such unclever algorithm usable.
const int N = 5;

vec2  bp[N]; // center of bubble
float br[N]; // radius of bubble
float bd[N]; // signed distance to bubble

float dbubbles(vec2 q, out int bid)
{
    float h = 3e38;
    bid = -1;
    // compute all the distances to query once to save some sqrts; can compute outside distance as we go, if negative, is inside at least one bubble
    for (int i = N + IZERO; i-- > 0; ) {
        vec2 a = bp[i]; 
        float r = br[i],
            ad = bd[i] = distance(q, a) - r;
        h = min(h, ad); // TODO smoothmin would make a meniscus
    }
    bool inside = h < 0.;
    if (inside)
    // for all bubbles
    for (int i = N + IZERO; i-- > 0; ) {
        float da = bd[i];
        // query is actually inside?
        if (da <= 0.)
        {
            // FIXME O(N^2), should at least go for O(log(N))
            // so far no luck with that, so doing it slowly
            for (int j = N; j-- > 0; ) {
                if (j == i) continue;
                float db = bd[j];
                // TODO optimization exclude query points on far side away from other bubble
                float bd = distance(bp[i], bp[j]); // FIXME still computed wastefully 2x each due to O(N^2) loop
                // bubbles touch?
                if (bd < br[i]+br[j]) {
                    // handle bubble completely inside another
                    if (bd < abs(br[i] - br[j])) continue;
                    float dq = dbubblesplit(q, bd, bp[i], bp[j], br[i], br[j]);
                    // if q in another bubble, give up on this one
                    if (dq < 0.) { da = 0.; break; }
                    da = max(da, -abs(dq));
                }
            }
            if (da < 0.) {
                bid = i;
                h = da;
                break;
            }
        }
    }
    return h;
}
// now separated out from the mainImage setup and drawing, ungolfed for clarity

void initBubbles()
{
    const float rn = 1./float(N);
    for (int i = N + IZERO; i-- > 0; ) {
        float z = float(i) * rn;
        bp[i] = vec2(0,0) + vec2(2.*sin(vec2(0,11) + .3*iTime + vec2(3,4) * .1*iTime + 6.*z));
        br[i] = 1.5 - z; //1.; //
    }
}

void mainImage(out vec4 o, vec2 p)
{
    vec2 r = iResolution.xy;
    float S = 8./r.y;
    vec2 q = (p - .5*r) * S;
    initBubbles();
    if (iMouse.z > 0.) bp[2] = (iMouse.xy - .5*r) * S;            
    // inside any bubble?
    int bid = -1;
    float h = dbubbles(q, bid);
    bool inside = h < 0.;
    h = abs(h);
    float T = 1.5*S; // half thickness
    h -= T;
    h /= S;
    float g = 1. - clamp(h, 0., 1.);
    vec3 c = inside ? vec3(1) : vec3(.1);
    const float phi = (sqrt(5.) + 1.) * .5;
    if (inside) // don't want to show outside voronoi regions here, didn't compute bid for those
        c = mix(c, cos(vec3(0,2,4) + phi * float(bid)) * .5 + .5, .5);
    c = mix(vec3(.2), c, mix(-cos(1.2*h), 1., .8)); // dim isocontours
    c = mix(c, vec3(1), g); // white outline
    o = vec4(pow(c, vec3(.45)), 1);
}

// uiop posted bubbles at http://shadertoy.com/view/DdlXWn before I could finish,
// which prompted me to public this way before it was even working properly.
// It's better now, but still needs much optimization and cleaning.

/*
    //r = r0*rd / (rd-r0); // TESTING
    //r = abs(r0*rd / (rd-r0)); // TESTING
    //r = 1./(1./r0 - 1./rd); // TESTING
    //r = 1./(1./rd - 1./r0); // TESTING
// so far seems doesn't work out the same,
// was hoping it would let me avoid the sqrt!
// TODO apparently 1/r2 = 1/r0-1/r1 or something like that,
// see question at mla's link where it's a given
// http://math.stackexchange.com/questions/680651/geometric-construction-of-2-soap-bubbles-meeting
// which seems to imply r2 = r0*r1 / (r1-r0)
// I think that may let me work backward to
// find the position of the smaller bubble that
// will produce 120 degree intersections.
// but doesn't appear to be a true relationship here.
// perhaps if actually moved the smaller bubble...
*/