Implicit Function for a Teapot

I haven't been solving many Project Euler problems for a long time, because it turns out that there is a lot of schoolwork to do when you take 5 courses, all CS or MATH, all 300-level or higher.

Anyway, something cool from one of my classes, an implicit function for a teapot, written in GLSL, with possibly a couple helper functions not included here. I may port this to a ShaderToy and/or improve the function at some point, but for now, here it is, and here is a video showing off the results:
https://www.youtube.com/watch?v=DSgwaYSnfcQ

float smin(float a, float b, float blendRadius) {
    float c = saturate(0.5 + (b - a) * (0.5 / blendRadius));
    return lerp(b, a, c) - blendRadius * c * (1.0 - c);
}

float teapotDistance(vec3 X) {

    //Spout - this is the most complex part
    vec3 C = vec3(-1.6f, -3.0f, 0);

    //Radius gets thinner until the end of the spout where it grows very quickly
    //fatness of the base of the spout is not dealt with here: it is increased
    //by using a large blend radius with the body
    float r = 0.7f +-.7 * sin(max(X.x, 3.0f));
    r = (X.x > 3.0f ? r + sin(X.x - 3.0f) : r);
    float e = 3.0f;

    //rotate a cylinder constantly to make a spout
    mat3 rot = rotation(pi / 2 + 2 * sin(X.x / 5.0f), 0, 0);
    vec3 tX = rot * X;
    tX.x -= .5f;
    vec2 d = abs(vec2(length(tX.xz - C.xz), tX.y - C.y)) - vec2(r, e);

    //Subtract a box in order to flatten the top of the spout
    C = vec3(2, 4.5f, 0);
    vec3 b = vec3(5.0f, 6.0f, 2.0f) - C;
    vec3 d3 = abs(X - C) - b;
    float box = min(maxComponent(d3), 0) + length(max(d3, vec3(0)));
    float spout = max(min(maxComponent(d), 0) + length(max(d, vec2(0, 0))), -box);

    //Body is another warped cylinder: it has no rotation, but its radius varies
    //as a square cosine (this results in a smaller curvature than the sine used in the spout)
    r = 3.0f + square(cos(X.y / 2.5f));
    e = 2.3f;
    C = vec3(-2, 1.0f, 0);
    d = abs(vec2(length(X.xz - C.xz), X.y - C.y)) - vec2(r, e);
    float body = min(maxComponent(d), 0) + length(max(d, vec2(0, 0)));


    //The bottom of the teapot is another warped cylinder, it is made to agree with
    //the bottom of the body, but then it shrinks rapidly before being clamped in order
    //to make the flat bottom
    float y = max(X.y, -1.8f);
    r = 3.0f + square(cos(-1.3f / 2.5f)) - square((y + 1.3f) / (.5f));
    e = .8f;
    C = vec3(-2, -1.3f, 0);
    d = abs(vec2(length(X.xz - C.xz), X.y - C.y)) - vec2(r, e);
    float bottom = min(maxComponent(d), 0) + length(max(d, vec2(0, 0)));

    //The top lid of the teapot is very similar to the bottom,
    //except that the clamp happens later to make the shaft of the top handle
    y = min(X.y, 3.81f);
    r = 3.0f + square(cos(3.3f / 2.5f)) - square((y - 3.3f)) / .1f;
    e = 0.9f;
    C = vec3(-2, 3.3f, 0);
    d = abs(vec2(length(X.xz - C.xz), X.y - C.y)) - vec2(r, e);
    float top = min(maxComponent(d), 0) + length(max(d, vec2(0, 0)));

    //The top of the handle of the lid is formed by making a hemisphere
    //A hemisphere is made by intersecting a sphere with a box
    C = vec3(-2, 4.4f, 0);
    r = 0.6f;
    float toptop = length(X - C) - r;
    C = vec3(-2, 3.3f, 0);
    b = vec3(-1.0f, 4.4f, 2.0f) - C;
    d3 = abs(X - C) - b;
    float boxtop = min(maxComponent(d3), 0) + length(max(d3, vec3(0)));
    toptop = max(toptop, boxtop);
    
    //The handle is a warped torus: the center is moved up along the length of the
    //torus, causing the bend
    float R = 0.2f;
    r = 1.9f;
    C = vec3(-2.0f + X.x/2.0f, 1.0f, 0.0f);
    tX = X;
    tX.y += (X.x + 7.0f)/3.0f - 0.5f;
    float handle = length(vec2(length(tX.xy - C.xy) - r, tX.z - C.z)) - R;
    //We blend everything together, most blend factors are just small epsilons
    //However, the body/spout blend has a large radius to fatten the bottom

    return smin(toptop, smin(top, smin(bottom, smin(handle, smin(spout, body, 0.5f), 0.001f), 0.001f), 0.001f), 0.001f); 
}