        const float J = bhAngularMomentum;
        const float M = bhMass;
        const float a = J / (M * c);

        mat4[4] grad;
        grad[0] = mat4(0.0);

        // `r` is defined by `((x^2 + y^2) / (r^2 + a^2)) + (z^2) / (r^2) = 1`
        // This came from: https://michaelmoroz.github.io/TracingGeodesics/
        float rho = dot(x.yzw, x.yzw) - a * a;
        float rho_dy = 2.0 * x.y;
        float rho_dz = 2.0 * x.z;
        float rho_dw = 2.0 * x.w;
        float r2_rootp_inner = rho * rho + 4.0 * a * a * x.w * x.w;
        float r2_rootp_inner_dy = 2.0 * rho_dy * rho;
        float r2_rootp_inner_dz = 2.0 * rho_dz * rho;
        float r2_rootp_inner_dw = 2.0 * rho_dw * rho + 8.0 * a * a * x.w;
        float r2_rootp = sqrt(r2_rootp_inner);
        float r2_rootp_dy = 0.5 * r2_rootp_inner_dy * inversesqrt(r2_rootp_inner);
        float r2_rootp_dz = 0.5 * r2_rootp_inner_dz * inversesqrt(r2_rootp_inner);
        float r2_rootp_dw = 0.5 * r2_rootp_inner_dw * inversesqrt(r2_rootp_inner);
        float r2 = 0.5 * (rho + r2_rootp);
        float r2_dy = 0.5 * (rho_dy + r2_rootp_dy);
        float r2_dz = 0.5 * (rho_dz + r2_rootp_dz);
        float r2_dw = 0.5 * (rho_dw + r2_rootp_dw);
        float r = sqrt(r2);
        float r_dy = 0.5 * r2_dy * inversesqrt(r2);
        float r_dz = 0.5 * r2_dz * inversesqrt(r2);
        float r_dw = 0.5 * r2_dw * inversesqrt(r2);

        float k_x_num = r * x.y + a * x.z;
        float k_x_num_dy = r_dy * x.y + r;
        float k_x_num_dz = r_dz * x.y + a;
        float k_x_num_dw = r_dw * x.y;
        float k_x_den = r2 + a * a;
        float k_x     = k_x_num / k_x_den;
        float k_x_dy  = (k_x_num_dy * k_x_den - k_x_num * r2_dy) / (k_x_den * k_x_den);
        float k_x_dz  = (k_x_num_dz * k_x_den - k_x_num * r2_dz) / (k_x_den * k_x_den);
        float k_x_dw  = (k_x_num_dw * k_x_den - k_x_num * r2_dw) / (k_x_den * k_x_den);
        if (isnan(k_x) || isinf(k_x)) { k_x = 0.0; }
        if (isnan(k_x_dy) || isinf(k_x_dy)) { k_x_dy = 0.0; }
        if (isnan(k_x_dz) || isinf(k_x_dz)) { k_x_dz = 0.0; }
        if (isnan(k_x_dw) || isinf(k_x_dw)) { k_x_dw = 0.0; }

        float k_y_num = r * x.z - a * x.y;
        float k_y_num_dy = r_dy * x.z - a;
        float k_y_num_dz = r_dz * x.z + r;
        float k_y_num_dw = r_dw * x.z;
        float k_y_den = r2 + a * a;
        float k_y     = k_y_num / k_y_den;
        float k_y_dy  = (k_y_num_dy * k_y_den - k_y_num * r2_dy) / (k_y_den * k_y_den);
        float k_y_dz  = (k_y_num_dz * k_y_den - k_y_num * r2_dz) / (k_y_den * k_y_den);
        float k_y_dw  = (k_y_num_dw * k_y_den - k_y_num * r2_dw) / (k_y_den * k_y_den);
        if (isnan(k_y) || isinf(k_y)) { k_y = 0.0; }
        if (isnan(k_y_dy) || isinf(k_y_dy)) { k_y_dy = 0.0; }
        if (isnan(k_y_dz) || isinf(k_y_dz)) { k_y_dz = 0.0; }
        if (isnan(k_y_dw) || isinf(k_y_dw)) { k_y_dw = 0.0; }

        float k_z     = x.w / r;
        float k_z_dy  = -x.w * r_dy / r2;
        float k_z_dz  = -x.w * r_dz / r2;
        float k_z_dw  = (r - x.w * r_dw) / r2;
        if (isnan(k_z) || isinf(k_z)) { k_z = 0.0; }
        if (isnan(k_z_dy) || isinf(k_z_dy)) { k_z_dy = 0.0; }
        if (isnan(k_z_dz) || isinf(k_z_dz)) { k_z_dz = 0.0; }
        if (isnan(k_z_dw) || isinf(k_z_dw)) { k_z_dw = 0.0; }

        vec3 k_vec    = vec3(k_x, k_y, k_z);
        vec3 k_vec_dy = vec3(k_x_dy, k_y_dy, k_z_dy);
        vec3 k_vec_dz = vec3(k_x_dz, k_y_dz, k_z_dz);
        vec3 k_vec_dw = vec3(k_x_dw, k_y_dw, k_z_dw);

        mat4 ku_kv = mat4(
            vec4(1.0, k_vec),
            vec4(k_vec.x, k_vec.x * k_vec),
            vec4(k_vec.y, k_vec.y * k_vec),
            vec4(k_vec.z, k_vec.z * k_vec)
        );
        mat4 ku_kv_dy = mat4(
            vec4(0.0, k_vec_dy),
            vec4(k_vec_dy.x, k_vec_dy.x * k_vec + k_vec.x * k_vec_dy),
            vec4(k_vec_dy.y, k_vec_dy.y * k_vec + k_vec.y * k_vec_dy),
            vec4(k_vec_dy.z, k_vec_dy.z * k_vec + k_vec.z * k_vec_dy)
        );
        mat4 ku_kv_dz = mat4(
            vec4(0.0, k_vec_dz),
            vec4(k_vec_dz.x, k_vec_dz.x * k_vec + k_vec.x * k_vec_dz),
            vec4(k_vec_dz.y, k_vec_dz.y * k_vec + k_vec.y * k_vec_dz),
            vec4(k_vec_dz.z, k_vec_dz.z * k_vec + k_vec.z * k_vec_dz)
        );
        mat4 ku_kv_dw = mat4(
            vec4(0.0, k_vec_dw),
            vec4(k_vec_dw.x, k_vec_dw.x * k_vec + k_vec.x * k_vec_dw),
            vec4(k_vec_dw.y, k_vec_dw.y * k_vec + k_vec.y * k_vec_dw),
            vec4(k_vec_dw.z, k_vec_dw.z * k_vec + k_vec.z * k_vec_dw)
        );

        float f_num = 2.0 * M * r2 * r - r2;
        float f_num_dy = 2.0 * M * (r2_dy * r + r2 * r_dy) - r2_dy;
        float f_num_dz = 2.0 * M * (r2_dz * r + r2 * r_dz) - r2_dz;
        float f_num_dw = 2.0 * M * (r2_dw * r + r2 * r_dw) - r2_dw;
        float f_den = r2 * r2 + a * a * x.w * x.w;
        float f_den_dy = 2.0 * r2_dy * r2;
        float f_den_dz = 2.0 * r2_dz * r2;
        float f_den_dw = 2.0 * r2_dw * r2 + 2.0 * a * a * x.w;
        float f = f_num / f_den;
        float f_dy = (f_num_dy * f_den - f_num * f_den_dy) / (f_den * f_den);
        float f_dz = (f_num_dz * f_den - f_num * f_den_dz) / (f_den * f_den);
        float f_dw = (f_num_dw * f_den - f_num * f_den_dw) / (f_den * f_den);

        if (isnan(f) || isinf(f)) { f = 0.0; }
        if (isnan(f_dy) || isinf(f_dy)) { f_dy = 0.0; }
        if (isnan(f_dz) || isinf(f_dz)) { f_dz = 0.0; }
        if (isnan(f_dw) || isinf(f_dw)) { f_dw = 0.0; }