Function ArcLength3

src/SB/Core/x/xSpline.cpp:161–218

Implementation of Composite Simpson's 1/3 Rule to calculate arc length

// Implementation of Composite Simpson's 1/3 Rule to calculate arc length
F32 ArcLength3(xCoef3* coef, F64 ustart, F64 uend)
{
    U32 i;
    F64 u;
    F64 h;
    F64 sum;

    F64 A;
    F64 B;
    F64 C;
    F64 D;
    F64 E;

    F64 u_eval;
    F64 temp_y1;
    F64 temp_z1;

    A = (coef->x).a[0];
    C = (coef->y).a[0];
    h = (coef->z).a[0];

    B = (coef->x).a[1];
    temp_y1 = (coef->y).a[1];
    temp_z1 = (coef->z).a[1];

    sum = (coef->x).a[2];
    u = (coef->y).a[2];
    u_eval = (coef->z).a[2];

    E = (h * h + A * A + C * C) * 9.0;
    D = (h * temp_z1 + A * B + C * temp_y1) * 12.0;

    C = (temp_z1 * temp_z1 + B * B + temp_y1 * temp_y1) * 4.0 +
        (h * u_eval + A * sum + C * u) * 6.0;
    A = u_eval * u_eval + sum * sum + u * u;
    h = (uend - ustart) / 50.0;
    B = (temp_z1 * u_eval + B * sum + temp_y1 * u) * 4.0;
    u = ustart + h;
    sum = 0.0;

    for (i = 2; i <= 51; i += 1)
    {
        if ((i & 1) == 0)
        {
            u_eval = xsqrt(A + u * (B + u * (C + u * (D + E * u)))) * 4.0;
        }
        else
        {
            u_eval = xsqrt(A + u * (B + u * (C + u * (D + E * u)))) * 2.0;
        }
        sum = sum + u_eval;
        u = u + h;
    }

    return (h * (sum + xsqrt(A + ustart * (B + ustart * (C + ustart * (D + E * ustart)))) +
                 xsqrt(A + uend * (B + uend * (C + uend * (D + E * uend)))))) /
           3.0;
}