Pen Settings

HTML

About HTML Preprocessors

HTML preprocessors can make writing HTML more powerful or convenient. For instance, Markdown is designed to be easier to write and read for text documents and you could write a loop in Pug.

Learn more · Versions

Adding Classes

In CodePen, whatever you write in the HTML editor is what goes within the <body> tags in a basic HTML5 template. So you don't have access to higher-up elements like the <html> tag. If you want to add classes there that can affect the whole document, this is the place to do it.

About the <head>

In CodePen, whatever you write in the HTML editor is what goes within the <body> tags in a basic HTML5 template. If you need things in the <head> of the document, put that code here.

!
Insecure Resource

The resource you are linking to is using the 'http' protocol, which may not work when the browser is using https.

CSS

About CSS Preprocessors

CSS preprocessors help make authoring CSS easier. All of them offer things like variables and mixins to provide convenient abstractions.

Learn more · Versions

CSS Base

About CSS Base

It's a common practice to apply CSS to a page that styles elements such that they are consistent across all browsers. We offer two of the most popular choices: normalize.css and a reset. Or, choose Neither and nothing will be applied.

Vendor Prefixing

About Vendor Prefixing

To get the best cross-browser support, it is a common practice to apply vendor prefixes to CSS properties and values that require them to work. For instance -webkit- or -moz-.

We offer two popular choices: Autoprefixer (which processes your CSS server-side) and -prefix-free (which applies prefixes via a script, client-side).

Add External Stylesheets/Pens

Any URLs added here will be added as <link>s in order, and before the CSS in the editor. You can use the CSS from another Pen by using its URL and the proper URL extension.

Powered by Algolia
About External Resources

You can apply CSS to your Pen from any stylesheet on the web. Just put a URL to it here and we'll apply it, in the order you have them, before the CSS in the Pen itself.

You can also link to another Pen here (use the .css URL Extension) and we'll pull the CSS from that Pen and include it. If it's using a matching preprocessor, use the appropriate URL Extension and we'll combine the code before preprocessing, so you can use the linked Pen as a true dependency.

Learn more

+ add another resource

JavaScript

About JavaScript Preprocessors

JavaScript preprocessors can help make authoring JavaScript easier and more convenient.

Learn more · Versions

Babel includes JSX processing.

Add External Scripts/Pens

Any URL's added here will be added as <script>s in order, and run before the JavaScript in the editor. You can use the URL of any other Pen and it will include the JavaScript from that Pen.

Powered by Algolia
About External Resources

You can apply a script from anywhere on the web to your Pen. Just put a URL to it here and we'll add it, in the order you have them, before the JavaScript in the Pen itself.

If the script you link to has the file extension of a preprocessor, we'll attempt to process it before applying.

You can also link to another Pen here, and we'll pull the JavaScript from that Pen and include it. If it's using a matching preprocessor, we'll combine the code before preprocessing, so you can use the linked Pen as a true dependency.

Learn more

+ add another resource

Packages

Add Packages

Search for and use JavaScript packages from npm here. By selecting a package, an import statement will be added to the top of the JavaScript editor for this package.

Powered by Algolia
About Packages

Using packages here is powered by Skypack, which makes packages from npm not only available on a CDN, but prepares them for native JavaScript ES6 import usage.

All packages are different, so refer to their docs for how they work.

If you're using React / ReactDOM, make sure to turn on Babel for the JSX processing.

Behavior

Save Automatically?

If active, Pens will autosave every 30 seconds after being saved once.

Auto-Updating Preview

If enabled, the preview panel updates automatically as you code. If disabled, use the "Run" button to update.

Format on Save

If enabled, your code will be formatted when you actively save your Pen. Note: your code becomes un-folded during formatting.

Editor Settings

Code Indentation

Want to change your Syntax Highlighting theme, Fonts and more?

Visit your global Editor Settings.

HTML 

              
                <script id="vertexShader" type="x-shader/x-vertex">

    void main(){
        gl_Position = projectionMatrix * viewMatrix * modelMatrix * vec4(position, 1.0);
    }

</script>

<script id="fragmentShader" type="x-shader/x-fragment">

    #define SDF_DERIVATIVE_TYPE 0
    #define NORMAL_DERIVATIVE_TYPE 0

    // ================ variables ================ //

    // ---------------- application ---------------- //

    struct App
    {
        float time;
        vec2 resolution;
    };

    struct Camera
    {
        vec3 position;
        mat4 viewMatrix;
        mat4 projectionMatrix;
    };

    struct Params
    {
        int numIterations;
        float convergenceCriteria;
        float finiteDifferenceEpsilon;
    };

    // ---------------- lighting ---------------- //

    struct PointLight
    {
        vec3 position;
        vec3 ambientColor;
        vec3 diffuseColor;
        vec3 specularColor;
    };

    struct DirectionalLight
    {
        vec3 direction;
        vec3 ambientColor;
        vec3 diffuseColor;
        vec3 specularColor;
    };

    struct Material
    {
        vec3 ambientColor;
        vec3 diffuseColor;
        vec3 specularColor;
        vec3 emissionColor;
        float shininess;
    };

    // ---------------- primitives ---------------- //

    struct Line
    {
        vec3 position;
        vec3 direction;
    };

    struct Sphere
    {
        vec3 position;
        float radius;
    };

    struct Intersection
    {
        vec3 position;
        bool intersected;
    };

    // ---------------- scene ---------------- //

    const int numLights = 2;

    struct Fractal
    {
        int power;
        int numIterations;
        float escapeCriteria;
    };

    struct Scene
    {
        vec3 backgroundColor;
        DirectionalLight lights[numLights];
        Material material;
        Sphere bound;
        Fractal fractal;
    };

    // ---------------- uniform ---------------- //

    uniform App uApp;
    uniform Camera uCamera;
    uniform Params uParams;
    uniform Scene uScene;

    // ================ functions ================ //

    // ---------------- utilities ---------------- //

    vec2 linmap(vec2 in_val, vec2 in_min, vec2 in_max, vec2 out_min, vec2 out_max)
    {
        return (in_val - in_min) / (in_max - in_min) * (out_max - out_min) + out_min;
    }

    // ---------------- primitives ---------------- //

    float sdfSphere(Sphere sphere, vec3 position)
    {
        return length(position - sphere.position) - sphere.radius;
    }

    vec3 normalSphere(Sphere sphere, vec3 position)
    {
        return normalize(position - sphere.position);
    }

    Intersection intersectionSphereLine(Sphere sphere, Line line)
    {
        vec3 difference = line.position - sphere.position;
        float a = dot(line.direction, line.direction);
        float b = 2.0 * dot(difference, line.direction);
        float c = dot(difference, difference) - pow(sphere.radius, 2.0);
        float d = pow(b, 2.0) - 4.0 * a * c;
        float t = (-b - sqrt(d)) / (2.0 * a);
        return Intersection(line.position + t * line.direction, d >= 0.0);
    }

    // ---------------- constructive solid geometry ---------------- //

    float csgUnion(float sd1, float sd2) { return min(sd1, sd2); }

    float csgSubtraction(float sd1, float sd2) { return max(-sd1, sd2); }

    float csgIntersection(float sd1, float sd2) { return max(sd1, sd2); }

    // ---------------- complex ---------------- //

    vec2 cAdd(vec2 c1, vec2 c2)
    {
        // return vec2(c1.x + c2.x, c1.y + c2.y);
        return c1 + c2;
    }

    vec2 cSub(vec2 c1, vec2 c2)
    {
        // return vec2(c1.x - c2.x, c1.y - c2.y);
        return c1 - c2;
    }

    vec2 cMul(vec2 c1, vec2 c2)
    {
        return vec2(c1.x * c2.x - c1.y * c2.y, c1.y * c2.x + c1.x * c2.y);
    }

    vec2 cConj(vec2 c)
    {
        return vec2(c.x, -c.y);
    }

    float cNorm(vec2 c)
    {
        // return sqrt(cMul(c, cConj(c)).x);
        return length(c);
    }

    vec2 cInv(vec2 c)
    {
        return cConj(c) / pow(cNorm(c), 2.0);
    }

    vec2 cDiv(vec2 c1, vec2 c2)
    {
        return cMul(c1, cInv(c2));
    }

    vec2 cPow(vec2 c, int n)
    {
        vec2 p = vec2(1.0, 0.0);
        for (int i = 0; i < n; ++i)
        {
            p = cMul(p, c);
        }
        return p;
    }

    // ---------------- quaternion ---------------- //

    vec4 qAdd(vec4 q1, vec4 q2)
    {
        // return vec4(q1.x + q2.x, q1.yzw + q2.yzw);
        return q1 + q2;
    }

    vec4 qSub(vec4 q1, vec4 q2)
    {
        // return vec4(q1.x - q2.x, q1.yzw - q2.yzw);
        return q1 - q2;
    }

    vec4 qMul(vec4 q1, vec4 q2)
    {
        return vec4(q1.x * q2.x - dot(q1.yzw, q2.yzw), q2.x * q1.yzw + q1.x * q2.yzw + cross(q1.yzw, q2.yzw));
    }

    vec4 qConj(vec4 q)
    {
        return vec4(q.x, -q.yzw);
    }

    float qNorm(vec4 q)
    {
        // return sqrt(qMul(q, qConj(q)).x);
        return length(q);
    }

    vec4 qInv(vec4 q)
    {
        return qConj(q) / pow(qNorm(q), 2.0);
    }

    vec4 qDiv(vec4 q1, vec4 q2)
    {
        return qMul(q1, qInv(q2));
    }

    vec4 qPow(vec4 q, int n)
    {
        vec4 p = vec4(1.0, vec3(0.0));
        for (int i = 0; i < n; ++i)
        {
            p = qMul(p, q);
        }
        return p;
    }

    // ---------------- trinion ---------------- //

    vec3 tAdd(vec3 t1, vec3 t2)
    {
        return t1 + t2;
    }

    vec3 tSub(vec3 t1, vec3 t2)
    {
        return t1 - t2;
    }

    vec3 tMul(vec3 t1, vec3 t2)
    {
        float r1 = length(t1);
        float r2 = length(t2);

        if (r1 > 0.0 && r2 > 0.0)
        {
            float a1 = asin(t1.z / r1);
            float a2 = asin(t2.z / r2);

            float b1 = atan(t1.y, t1.x);
            float b2 = atan(t2.y, t2.x);

            float r = r1 * r2;
            float a = a1 + a2;
            float b = b1 + b2;

            float x = r * cos(a) * cos(b);
            float y = r * cos(a) * sin(b);
            float z = r * sin(a);

            return vec3(x, y, z);
        }
        else
        {
            return vec3(0.0);
        }
    }

    vec3 tPow(vec3 t, int n)
    {
        /* NOTE: This does not work

        vec3 p = vec3(1.0, vec2(0.0));
        for (int i = 0; i < n; ++i)
        {
            p = tMul(p, t);
        }
        return p;

        */

        float r = length(t);

        if (r > 0.0)
        {
            float a = asin(t.z / r);
            float b = atan(t.y, t.x);

            float pr = pow(r, float(n));
            float pa = a * float(n);
            float pb = b * float(n);

            float x = pr * cos(pa) * cos(pb);
            float y = pr * cos(pa) * sin(pb);
            float z = pr * sin(pa);

            return vec3(x, y, z);
        }
        else
        {
            return vec3(0.0);
        }
    }

    // ---------------- dual ---------------- //

    struct DualQ
    {
        vec4 q;
        vec4 d;
    };

    DualQ dqAdd(DualQ dq1, DualQ dq2)
    {
        return DualQ(qAdd(dq1.q, dq2.q), qAdd(dq1.d, dq2.d));
    }

    DualQ dqSub(DualQ dq1, DualQ dq2)
    {
        return DualQ(qSub(dq1.q, dq2.q), qSub(dq1.d, dq2.d));
    }

    DualQ dqMul(DualQ dq1, DualQ dq2)
    {
        return DualQ(qMul(dq1.q, dq2.q), qAdd(qMul(dq1.d, dq2.q), qMul(dq1.q, dq2.d)));
    }

    DualQ dqDiv(DualQ dq1, DualQ dq2)
    {
        return DualQ(qDiv(dq1.q, dq2.q), qDiv(qSub(qMul(dq1.d, dq2.q), qMul(dq1.q, dq2.d)), qMul(dq2.q, dq2.q)));
    }

    DualQ dqPow(DualQ dq, int n)
    {
        DualQ dp = DualQ(vec4(1.0, vec3(0.0)), vec4(0.0, vec3(0.0)));
        for (int i = 0; i < n; ++i)
        {
            dp = dqMul(dp, dq);
        }
        return dp;
    }

    struct DualS
    {
        float s;
        float d;
    };

    DualS dAdd(DualS d1, DualS d2)
    {
        return DualS(d1.s + d2.s, d1.d + d2.d);
    }

    DualS dSub(DualS d1, DualS d2)
    {
        return DualS(d1.s - d2.s, d1.d - d2.d);
    }

    DualS dMul(DualS d1, DualS d2)
    {
        return DualS(d1.s * d2.s, d1.d * d2.s + d1.s * d2.d);
    }

    DualS dDiv(DualS d1, DualS d2)
    {
        return DualS(d1.s / d2.s, (d1.d * d2.s - d1.s * d2.d) / (d2.s * d2.s));
    }

    DualS dPow(DualS d, float n)
    {
        return DualS(pow(d.s, n), d.d * n * pow(d.s, n - 1.0));
    }

    DualS dSqrt(DualS d)
    {
        return DualS(sqrt(d.s), d.d * 0.5 / sqrt(d.s));
    }

    DualS dSin(DualS d)
    {
        return DualS(sin(d.s), d.d * cos(d.s));
    }

    DualS dCos(DualS d)
    {
        return DualS(cos(d.s), d.d * -sin(d.s));
    }

    DualS dTan(DualS d)
    {
        return DualS(tan(d.s), d.d / pow(cos(d.s), 2.0));
    }

    DualS dArcSin(DualS d)
    {
        return DualS(asin(d.s), d.d / sqrt(1.0 - pow(d.s, 2.0)));
    }

    DualS dArcCos(DualS d)
    {
        return DualS(acos(d.s), d.d / -sqrt(1.0 - pow(d.s, 2.0)));
    }

    DualS dArcTan(DualS d)
    {
        return DualS(atan(d.s), d.d / (1.0 + pow(d.s, 2.0)));
    }

    struct DualT
    {
        vec3 t;
        vec3 d;
    };

    DualT dtAdd(DualT dt1, DualT dt2)
    {
        return DualT(dt1.t + dt2.t, dt1.d + dt2.d);
    }

    DualT dtSub(DualT dt1, DualT dt2)
    {
        return DualT(dt1.t - dt2.t, dt1.d - dt2.d);
    }

    DualT dtMul(DualT dt1, DualT dt2)
    {
        // return DualT(tMul(dt1.t, dt2.t), tAdd(tMul(dt1.d, dt2.t), tMul(dt1.t, dt2.d)));

        DualS dx1 = DualS(dt1.t.x, dt1.d.x);
        DualS dy1 = DualS(dt1.t.y, dt1.d.y);
        DualS dz1 = DualS(dt1.t.z, dt1.d.z);

        DualS dx2 = DualS(dt2.t.x, dt2.d.x);
        DualS dy2 = DualS(dt2.t.y, dt2.d.y);
        DualS dz2 = DualS(dt2.t.z, dt2.d.z);

        DualS dr1 = dSqrt(dAdd(dPow(dx1, 2.0), dAdd(dPow(dy1, 2.0), dPow(dz1, 2.0))));
        DualS dr2 = dSqrt(dAdd(dPow(dx2, 2.0), dAdd(dPow(dy2, 2.0), dPow(dz2, 2.0))));

        if (dr1.s > 0.0 && dr2.s > 0.0)
        {
            DualS da1 = dArcSin(dDiv(dz1, dr1));
            DualS da2 = dArcSin(dDiv(dz2, dr2));

            DualS db1 = dArcTan(dDiv(dy1, dx1));
            DualS db2 = dArcTan(dDiv(dy2, dx2));

            DualS dr = dMul(dr1, dr2);
            DualS da = dAdd(da1, da2);
            DualS db = dAdd(db1, db2);

            DualS dx = dMul(dr, dMul(dCos(da), dCos(db)));
            DualS dy = dMul(dr, dMul(dCos(da), dSin(db)));
            DualS dz = dMul(dr, dSin(da));

            return DualT(vec3(dx.s, dy.s, dz.s), vec3(dx.d, dy.d, dz.d));
        }
        else
        {
            return DualT(vec3(0.0), vec3(0.0));
        }
    }

    DualT dtPow(DualT dt, int n)
    {
        /* NOTE: This does not work

        DualT dp = DualT(vec3(1.0, vec2(0.0)), vec3(0.0, vec2(0.0)));
        for (int i = 0; i < n; ++i)
        {
            dp = dtMul(dp, dt);
        }
        return dp;

        */

        DualS dx = DualS(dt.t.x, dt.d.x);
        DualS dy = DualS(dt.t.y, dt.d.y);
        DualS dz = DualS(dt.t.z, dt.d.z);

        DualS dr = dSqrt(dAdd(dPow(dx, 2.0), dAdd(dPow(dy, 2.0), dPow(dz, 2.0))));

        if (dr.s > 0.0)
        {
            DualS da = dArcSin(dDiv(dz, dr));
            DualS db = dArcTan(dDiv(dy, dx));

            DualS dpr = dPow(dr, float(n));
            DualS dpa = dMul(da, DualS(float(n), 0.0));
            DualS dpb = dMul(db, DualS(float(n), 0.0));

            DualS dx = dMul(dpr, dMul(dCos(dpa), dCos(dpb)));
            DualS dy = dMul(dpr, dMul(dCos(dpa), dSin(dpb)));
            DualS dz = dMul(dpr, dSin(dpa));

            return DualT(vec3(dx.s, dy.s, dz.s), vec3(dx.d, dy.d, dz.d));
        }
        else
        {
            return DualT(vec3(0.0), vec3(0.0));
        }
    }

    // ---------------- fractals ---------------- //

    #if SDF_DERIVATIVE_TYPE == 0
    float sdfJulia(Fractal fractal, vec4 z, vec4 c)
    {
        DualQ dzx = DualQ(z, vec4(1.0, 0.0, 0.0, 0.0));
        DualQ dzy = DualQ(z, vec4(0.0, 1.0, 0.0, 0.0));
        DualQ dzz = DualQ(z, vec4(0.0, 0.0, 1.0, 0.0));
        DualQ dzw = DualQ(z, vec4(0.0, 0.0, 0.0, 1.0));
        
        DualQ dcx = DualQ(c, vec4(0.0, 0.0, 0.0, 0.0));
        DualQ dcy = DualQ(c, vec4(0.0, 0.0, 0.0, 0.0));
        DualQ dcz = DualQ(c, vec4(0.0, 0.0, 0.0, 0.0));
        DualQ dcw = DualQ(c, vec4(0.0, 0.0, 0.0, 0.0));
        
        for (int i = 0; i < fractal.numIterations; ++i)
        {
            // forward-mode automatic differentiation
            dzx = dqAdd(dqPow(dzx, fractal.power), dcx);
            dzy = dqAdd(dqPow(dzy, fractal.power), dcy);
            dzz = dqAdd(dqPow(dzz, fractal.power), dcz);
            dzw = dqAdd(dqPow(dzw, fractal.power), dcw);
            
            if (qNorm(dzx.q) > fractal.escapeCriteria) break;
        }
        
        mat4 J = mat4(dzx.d, dzy.d, dzz.d, dzw.d);
        return (qNorm(dzx.q) * log(qNorm(dzx.q))) / (2.0 * qNorm(normalize(dzx.q) * J));
    }
    float sdfMandelbrot(Fractal fractal, vec4 c, vec4 z)
    {
        DualQ dzx = DualQ(z, vec4(0.0, 0.0, 0.0, 0.0));
        DualQ dzy = DualQ(z, vec4(0.0, 0.0, 0.0, 0.0));
        DualQ dzz = DualQ(z, vec4(0.0, 0.0, 0.0, 0.0));
        DualQ dzw = DualQ(z, vec4(0.0, 0.0, 0.0, 0.0));
        
        DualQ dcx = DualQ(c, vec4(1.0, 0.0, 0.0, 0.0));
        DualQ dcy = DualQ(c, vec4(0.0, 1.0, 0.0, 0.0));
        DualQ dcz = DualQ(c, vec4(0.0, 0.0, 1.0, 0.0));
        DualQ dcw = DualQ(c, vec4(0.0, 0.0, 0.0, 1.0));
        
        for (int i = 0; i < fractal.numIterations; ++i)
        {
            // forward-mode automatic differentiation
            dzx = dqAdd(dqPow(dzx, fractal.power), dcx);
            dzy = dqAdd(dqPow(dzy, fractal.power), dcy);
            dzz = dqAdd(dqPow(dzz, fractal.power), dcz);
            dzw = dqAdd(dqPow(dzw, fractal.power), dcw);
            
            if (qNorm(dzx.q) > fractal.escapeCriteria) break;
        }
        
        mat4 J = mat4(dzx.d, dzy.d, dzz.d, dzw.d);
        return (qNorm(dzx.q) * log(qNorm(dzx.q))) / (2.0 * qNorm(normalize(dzx.q) * J));
    }
    float sdfMandelbulb(Fractal fractal, vec3 c, vec3 z)
    {
        DualT dzx = DualT(z, vec3(0.0, 0.0, 0.0));
        DualT dzy = DualT(z, vec3(0.0, 0.0, 0.0));
        DualT dzz = DualT(z, vec3(0.0, 0.0, 0.0));
        
        DualT dcx = DualT(c, vec3(1.0, 0.0, 0.0));
        DualT dcy = DualT(c, vec3(0.0, 1.0, 0.0));
        DualT dcz = DualT(c, vec3(0.0, 0.0, 1.0));
        
        for (int i = 0; i < fractal.numIterations; ++i)
        {
            // forward-mode automatic differentiation
            dzx = dtAdd(dtPow(dzx, fractal.power), dcx);
            dzy = dtAdd(dtPow(dzy, fractal.power), dcy);
            dzz = dtAdd(dtPow(dzz, fractal.power), dcz);
            
            if (length(dzx.t) > fractal.escapeCriteria) break;
        }

        mat3 J = mat3(dzx.d, dzy.d, dzz.d);
        return (length(dzx.t) * log(length(dzx.t))) / (2.0 * length(normalize(dzx.t) * J));
    }
    #elif SDF_DERIVATIVE_TYPE == 1
    float sdfJulia(Fractal fractal, vec4 z, vec4 c)
    {
        vec4 dzx = vec4(1.0, 0.0, 0.0, 0.0);
        vec4 dzy = vec4(0.0, 1.0, 0.0, 0.0);
        vec4 dzz = vec4(0.0, 0.0, 1.0, 0.0);
        vec4 dzw = vec4(0.0, 0.0, 0.0, 1.0);
        
        vec4 dcx = vec4(0.0, 0.0, 0.0, 0.0);
        vec4 dcy = vec4(0.0, 0.0, 0.0, 0.0);
        vec4 dcz = vec4(0.0, 0.0, 0.0, 0.0);
        vec4 dcw = vec4(0.0, 0.0, 0.0, 0.0);
        
        for (int i = 0; i < fractal.numIterations; ++i)
        {
            vec4 zp = qPow(z, fractal.power - 1);
            
            // forward-mode manual differentiation
            dzx = qAdd(float(fractal.power) * qMul(zp, dzx), dcx);
            dzy = qAdd(float(fractal.power) * qMul(zp, dzy), dcy);
            dzz = qAdd(float(fractal.power) * qMul(zp, dzz), dcz);
            dzw = qAdd(float(fractal.power) * qMul(zp, dzw), dcw);
            
            z = qAdd(qMul(zp, z), c);
            
            if (qNorm(z) > fractal.escapeCriteria) break;
        }
        
        mat4 J = mat4(dzx, dzy, dzz, dzw);
        return (qNorm(z) * log(qNorm(z))) / (2.0 * qNorm(normalize(z) * J));
    }
    float sdfMandelbrot(Fractal fractal, vec4 c, vec4 z)
    {
        vec4 dzx = vec4(0.0, 0.0, 0.0, 0.0);
        vec4 dzy = vec4(0.0, 0.0, 0.0, 0.0);
        vec4 dzz = vec4(0.0, 0.0, 0.0, 0.0);
        vec4 dzw = vec4(0.0, 0.0, 0.0, 0.0);
        
        vec4 dcx = vec4(1.0, 0.0, 0.0, 0.0);
        vec4 dcy = vec4(0.0, 1.0, 0.0, 0.0);
        vec4 dcz = vec4(0.0, 0.0, 1.0, 0.0);
        vec4 dcw = vec4(0.0, 0.0, 0.0, 1.0);
        
        for (int i = 0; i < fractal.numIterations; ++i)
        {
            vec4 zp = qPow(z, fractal.power - 1);
            
            // forward-mode manual differentiation
            dzx = qAdd(float(fractal.power) * qMul(zp, dzx), dcx);
            dzy = qAdd(float(fractal.power) * qMul(zp, dzy), dcy);
            dzz = qAdd(float(fractal.power) * qMul(zp, dzz), dcz);
            dzw = qAdd(float(fractal.power) * qMul(zp, dzw), dcw);
            
            z = qAdd(qMul(zp, z), c);
            
            if (qNorm(z) > fractal.escapeCriteria) break;
        }
        
        mat4 J = mat4(dzx, dzy, dzz, dzw);
        return (qNorm(z) * log(qNorm(z))) / (2.0 * qNorm(normalize(z) * J));
    }
    #else
    #endif

    #if NORMAL_DERIVATIVE_TYPE == 0
    vec4 normalJulia(Fractal fractal, vec4 z, vec4 c)
    {
        DualQ dzx = DualQ(z, vec4(1.0, 0.0, 0.0, 0.0));
        DualQ dzy = DualQ(z, vec4(0.0, 1.0, 0.0, 0.0));
        DualQ dzz = DualQ(z, vec4(0.0, 0.0, 1.0, 0.0));
        DualQ dzw = DualQ(z, vec4(0.0, 0.0, 0.0, 1.0));
        
        DualQ dcx = DualQ(c, vec4(0.0, 0.0, 0.0, 0.0));
        DualQ dcy = DualQ(c, vec4(0.0, 0.0, 0.0, 0.0));
        DualQ dcz = DualQ(c, vec4(0.0, 0.0, 0.0, 0.0));
        DualQ dcw = DualQ(c, vec4(0.0, 0.0, 0.0, 0.0));
        
        for (int i = 0; i < fractal.numIterations; ++i)
        {
            // forward-mode automatic differentiation
            dzx = dqAdd(dqPow(dzx, fractal.power), dcx);
            dzy = dqAdd(dqPow(dzy, fractal.power), dcy);
            dzz = dqAdd(dqPow(dzz, fractal.power), dcz);
            dzw = dqAdd(dqPow(dzw, fractal.power), dcw);
            
            if (qNorm(dzx.q) > fractal.escapeCriteria) break;
        }
        
        mat4 J = mat4(dzx.d, dzy.d, dzz.d, dzw.d);
        return dzx.q * J;
    }
    vec4 normalMandelbrot(Fractal fractal, vec4 c, vec4 z)
    {
        DualQ dzx = DualQ(z, vec4(0.0, 0.0, 0.0, 0.0));
        DualQ dzy = DualQ(z, vec4(0.0, 0.0, 0.0, 0.0));
        DualQ dzz = DualQ(z, vec4(0.0, 0.0, 0.0, 0.0));
        DualQ dzw = DualQ(z, vec4(0.0, 0.0, 0.0, 0.0));
        
        DualQ dcx = DualQ(c, vec4(1.0, 0.0, 0.0, 0.0));
        DualQ dcy = DualQ(c, vec4(0.0, 1.0, 0.0, 0.0));
        DualQ dcz = DualQ(c, vec4(0.0, 0.0, 1.0, 0.0));
        DualQ dcw = DualQ(c, vec4(0.0, 0.0, 0.0, 1.0));
        
        for (int i = 0; i < fractal.numIterations; ++i)
        {
            // forward-mode automatic differentiation
            dzx = dqAdd(dqPow(dzx, fractal.power), dcx);
            dzy = dqAdd(dqPow(dzy, fractal.power), dcy);
            dzz = dqAdd(dqPow(dzz, fractal.power), dcz);
            dzw = dqAdd(dqPow(dzw, fractal.power), dcw);
            
            if (qNorm(dzx.q) > fractal.escapeCriteria) break;
        }
        
        mat4 J = mat4(dzx.d, dzy.d, dzz.d, dzw.d);
        return dzx.q * J;
    }
    vec3 normalMandelbulb(Fractal fractal, vec3 c, vec3 z)
    {
        DualT dzx = DualT(z, vec3(0.0, 0.0, 0.0));
        DualT dzy = DualT(z, vec3(0.0, 0.0, 0.0));
        DualT dzz = DualT(z, vec3(0.0, 0.0, 0.0));
        
        DualT dcx = DualT(c, vec3(1.0, 0.0, 0.0));
        DualT dcy = DualT(c, vec3(0.0, 1.0, 0.0));
        DualT dcz = DualT(c, vec3(0.0, 0.0, 1.0));
        
        for (int i = 0; i < fractal.numIterations; ++i)
        {
            // forward-mode automatic differentiation
            dzx = dtAdd(dtPow(dzx, fractal.power), dcx);
            dzy = dtAdd(dtPow(dzy, fractal.power), dcy);
            dzz = dtAdd(dtPow(dzz, fractal.power), dcz);
            
            if (length(dzx.t) > fractal.escapeCriteria) break;
        }
        
        mat3 J = mat3(dzx.d, dzy.d, dzz.d);
        return dzx.t * J;
    }
    #elif NORMAL_DERIVATIVE_TYPE == 1
    vec4 normalJulia(Fractal fractal, vec4 z, vec4 c)
    {
        vec4 dzx = vec4(1.0, 0.0, 0.0, 0.0);
        vec4 dzy = vec4(0.0, 1.0, 0.0, 0.0);
        vec4 dzz = vec4(0.0, 0.0, 1.0, 0.0);
        vec4 dzw = vec4(0.0, 0.0, 0.0, 1.0);
        
        vec4 dcx = vec4(0.0, 0.0, 0.0, 0.0);
        vec4 dcy = vec4(0.0, 0.0, 0.0, 0.0);
        vec4 dcz = vec4(0.0, 0.0, 0.0, 0.0);
        vec4 dcw = vec4(0.0, 0.0, 0.0, 0.0);
        
        for (int i = 0; i < fractal.numIterations; ++i)
        {
            vec4 zp = qPow(z, fractal.power - 1);
            
            // forward-mode manual differentiation
            dzx = qAdd(float(fractal.power) * qMul(zp, dzx), dcx);
            dzy = qAdd(float(fractal.power) * qMul(zp, dzy), dcy);
            dzz = qAdd(float(fractal.power) * qMul(zp, dzz), dcz);
            dzw = qAdd(float(fractal.power) * qMul(zp, dzw), dcw);
            
            z = qAdd(qMul(zp, z), c);
            
            if (qNorm(z) > fractal.escapeCriteria) break;
        }
        
        mat4 J = mat4(dzx, dzy, dzz, dzw);
        return z * J;
    }
    vec4 normalMandelbrot(Fractal fractal, vec4 c, vec4 z)
    {
        vec4 dzx = vec4(0.0, 0.0, 0.0, 0.0);
        vec4 dzy = vec4(0.0, 0.0, 0.0, 0.0);
        vec4 dzz = vec4(0.0, 0.0, 0.0, 0.0);
        vec4 dzw = vec4(0.0, 0.0, 0.0, 0.0);
        
        vec4 dcx = vec4(1.0, 0.0, 0.0, 0.0);
        vec4 dcy = vec4(0.0, 1.0, 0.0, 0.0);
        vec4 dcz = vec4(0.0, 0.0, 1.0, 0.0);
        vec4 dcw = vec4(0.0, 0.0, 0.0, 1.0);
        
        for (int i = 0; i < fractal.numIterations; ++i)
        {
            vec4 zp = qPow(z, fractal.power - 1);
            
            // forward-mode manual differentiation
            dzx = qAdd(float(fractal.power) * qMul(zp, dzx), dcx);
            dzy = qAdd(float(fractal.power) * qMul(zp, dzy), dcy);
            dzz = qAdd(float(fractal.power) * qMul(zp, dzz), dcz);
            dzw = qAdd(float(fractal.power) * qMul(zp, dzw), dcw);
            
            z = qAdd(qMul(zp, z), c);
            
            if (qNorm(z) > fractal.escapeCriteria) break;
        }
        
        mat4 J = mat4(dzx, dzy, dzz, dzw);
        return z * J;
    }
    #else
    #endif

    // ---------------- reflection ---------------- //

    vec3 phongReflection(vec3 surfaceNormal, vec3 viewDirection, DirectionalLight lights[numLights], Material material)
    {
        surfaceNormal = normalize(surfaceNormal);
        viewDirection = normalize(viewDirection);
        
        vec3 ambientColor = vec3(0.0);
        vec3 diffuseColor = vec3(0.0);
        vec3 specularColor = vec3(0.0);
        
        for (int i = 0; i < lights.length(); ++i)
        {
            vec3 lightDirection = normalize(lights[i].direction);
            vec3 reflectedDirection = reflect(lightDirection, surfaceNormal);
            float diffuseCoefficient = max(dot(-lightDirection, surfaceNormal), 0.0);
            float specularCoefficient = pow(max(dot(reflectedDirection, -viewDirection), 0.0), material.shininess);
            ambientColor += lights[i].ambientColor * material.ambientColor;
            diffuseColor += lights[i].diffuseColor * material.diffuseColor * diffuseCoefficient;
            specularColor += lights[i].specularColor * material.specularColor * specularCoefficient;
        }
        
        vec3 color = clamp(ambientColor + diffuseColor + specularColor + material.emissionColor, 0.0, 1.0);
        return color;
    }

    // ---------------- sphere tracing ---------------- //

    vec3 sphereTracing(App app, Scene scene, Params params, Line ray)
    {
        Intersection intersection = intersectionSphereLine(scene.bound, ray);
        
        if (intersection.intersected)
        {
            ray.position = intersection.position;
            
            // The hyperparameter from Inigo Quilez (https://www.shadertoy.com/view/MsfGRr)
            vec4 juliaType = 0.45 * cos(vec4(0.5, 3.9, 1.4, 1.1) + app.time * 0.15 * vec4(1.2, 1.7, 1.3, 2.5)) - vec4(0.3, 0.0, 0.0, 0.0);
            vec4 criticalPoint = vec4(0.0);
            
            for (int i = 0; i < params.numIterations; ++i)
            {
                float sd = sdfMandelbulb(scene.fractal, ray.position, criticalPoint.xyz);
                
                // ray marching
                ray.position += sd * ray.direction;
                
                // collision detection
                if (abs(sd) < params.convergenceCriteria)
                {
                    vec3 surfaceNormal = normalize(normalMandelbulb(scene.fractal, ray.position, criticalPoint.xyz));
                    
                    vec3 fragColor = phongReflection(surfaceNormal, ray.direction, scene.lights, scene.material);
                    
                    return fragColor;
                }
                
                if (sdfSphere(scene.bound, ray.position) > 0.0) break;
            }
        }
        
        return scene.backgroundColor;
    }

    // ---------------- main ---------------- //

    void main()
    {
        vec2 fragCoord = linmap(gl_FragCoord.xy, vec2(0, 0), uApp.resolution, vec2(-1.0, -1.0), vec2(1.0, 1.0));
        vec3 rayDirection = normalize(inverse(mat3(uCamera.projectionMatrix) * mat3(uCamera.viewMatrix)) * vec3(fragCoord, 1.0));
        
        Line ray = Line(uCamera.position, rayDirection);
        vec3 fragColor = sphereTracing(uApp, uScene, uParams, ray);
        gl_FragColor = vec4(fragColor, 1.0);
    }

</script>
!

CSS
!

JS 

              
                const init = () => {

    const renderer = new THREE.WebGLRenderer({ antialias:true });
    document.body.appendChild(renderer.domElement);

    const scene = new THREE.Scene();

    const orthographiCamera = new THREE.OrthographicCamera(window.innerWidth / -2.0, window.innerWidth / +2.0, window.innerHeight / +2.0, window.innerHeight / -2.0, 0.0, 1.0);
    const perspectiveCamera = new THREE.PerspectiveCamera(45.0, window.innerWidth / window.innerHeight, 0.1, 1000.0);
    const controls = new THREE.OrbitControls(perspectiveCamera, renderer.domElement);
    const clock = new THREE.Clock();

    perspectiveCamera.position.set(0.0, 0.0, 5.0);
    perspectiveCamera.lookAt(new THREE.Vector3(0.0, 0.0, 0.0));

    const geometry = new THREE.PlaneBufferGeometry(window.innerWidth, window.innerHeight);

    const uniforms = {
        uApp: {
            value: {
                time: clock.getElapsedTime(),
                resolution: new THREE.Vector2(window.innerWidth, window.innerHeight)
            }
        },
        uCamera: {
            value: {
                position: perspectiveCamera.position,
                viewMatrix: perspectiveCamera.matrixWorldInverse,
                projectionMatrix: perspectiveCamera.projectionMatrix
            }
        },
        uParams: {
            value: {
                numIterations: 300,
                convergenceCriteria: 0.0001,
                finiteDifferenceEpsilon: 0.0001
            }
        },
        uScene: {
            value: {
                backgroundColor: new THREE.Vector3(0.0, 0.0, 0.0),
                lights: [
                    {
                        direction: new THREE.Vector3(1.0, 1.0, 1.0),
                        ambientColor: new THREE.Vector3(1.0, 1.0, 1.0),
                        diffuseColor: new THREE.Vector3(1.0, 1.0, 1.0),
                        specularColor: new THREE.Vector3(1.0, 1.0, 1.0)
                    },
                    {
                        direction: new THREE.Vector3(-1.0, -1.0, -1.0),
                        ambientColor: new THREE.Vector3(1.0, 1.0, 1.0),
                        diffuseColor: new THREE.Vector3(1.0, 1.0, 1.0),
                        specularColor: new THREE.Vector3(1.0, 1.0, 1.0)
                    }
                ],
                material: {
                    ambientColor: new THREE.Vector3(0.05, 0.05, 0.05),
                    diffuseColor: new THREE.Vector3(0.5, 0.5, 0.5),
                    specularColor: new THREE.Vector3(1.0, 1.0, 1.0),
                    emissionColor: new THREE.Vector3(0.0, 0.0, 0.0),
                    shininess: 64.0
                },
                bound: {
                    position: new THREE.Vector3(0.0, 0.0, 0.0),
                    radius: 2.0
                },
                fractal: {
                    power: 8,
                    numIterations: 4,
                    escapeCriteria: 2.0
                }
            }
        }
    }

    const material = new THREE.ShaderMaterial({
        vertexShader: document.getElementById('vertexShader').textContent,
        fragmentShader: document.getElementById('fragmentShader').textContent,
        uniforms: uniforms
    });

    scene.add(new THREE.Mesh(geometry, material));

    const onWindowResize = (event) => {
        uniforms.uApp.value.resolution.x = window.innerWidth * window.devicePixelRatio;
        uniforms.uApp.value.resolution.y = window.innerHeight * window.devicePixelRatio;
        // NOTE: https://ics.media/tutorial-three/renderer_resize/
        renderer.setPixelRatio(window.devicePixelRatio);
        renderer.setSize(window.innerWidth, window.innerHeight);
        perspectiveCamera.aspect = window.innerWidth / window.innerHeight;
        perspectiveCamera.updateProjectionMatrix();
    }
    onWindowResize();
    window.addEventListener('resize', onWindowResize, false);

    const animate = () => {

        requestAnimationFrame(animate);

        const update = () => {
            controls.update();
            perspectiveCamera.lookAt(new THREE.Vector3(0.0, 0.0, 0.0));
            uniforms.uApp.value.time = clock.getElapsedTime();
            uniforms.uCamera.value.position = perspectiveCamera.position;
            uniforms.uCamera.value.viewMatrix = perspectiveCamera.matrixWorldInverse;
            uniforms.uCamera.value.projectionMatrix = perspectiveCamera.projectionMatrix;
        }

        update();

        renderer.render(scene, orthographiCamera);
    };

    animate();
}

window.addEventListener("load", init);
!
999px

Console