Translated lens undistortion algorithm giving wrong results at corners

Question

I am trying to adapt the lens undistortion from Oliver Kreylos' Kinect library to Unity, but at the corners the vertices end up in the wrong positions for any radial distortion coefficients below about -0.15, give or take depending on the coefficient. The 2nd (k0), 4th (k1), and 6th (k2) order radial distortion coefficients are 0.091, -0.269, 0.094 respectively for my Kinect.

Here is the resulting mesh with k0 set to -0.18 and the rest to 0

Here is the undistortion function translated to C#

public Vector2 Undistort(Vector2 distorted)
{
    /*******************************************************************
    Must run two-dimensional Newton-Raphson iteration to invert the full
    distortion correction formula.
    *******************************************************************/

    /* Start iteration at the distorted point: */
    Vector2 p = distorted;
    for (int iteration = 0; iteration < undistortMaxSteps; ++iteration)
    {
        /* Calculate the function value at p: */
        Vector2 d = p - center;
        float r2 = d.sqrMagnitude;
        float div = 1 + (radialDistortion[0] + (radialDistortion[1] + radialDistortion[2] * r2) * r2) * r2;
        Vector2 fp = new Vector2(
            center[0] + d[0] * div + 2 * tangentialDistortion[0] * d[0] * d[1] + tangentialDistortion[1] * (r2 + 2 * d[0] * d[0]) - distorted[0],
            center[1] + d[1] * div + tangentialDistortion[0] * (r2 + 2 * d[1] * d[1]) + 2 * tangentialDistortion[1] * d[0] * d[1] - distorted[1]
            );

        /* Bail out if close enough: */
        if (fp.sqrMagnitude < undistortMaxError)
            break;

        /* Calculate the function derivative at p: */
        float[,] fpd = new float[2, 2];
        fpd[0, 0] = div + d[0] * (2 * radialDistortion[0] + (4 * radialDistortion[1] + 6 * radialDistortion[2] * r2) * r2) * d[0] + 2 * tangentialDistortion[0] * d[1] + 6 * tangentialDistortion[1] * d[0]; // d fp[0] / d p[0]
        fpd[0, 1] = d[0] * (2 * radialDistortion[0] + (4 * radialDistortion[1] + 6 * radialDistortion[2] * r2) * r2) * d[1] + 2 * tangentialDistortion[0] * d[0] + 2 * tangentialDistortion[1] * d[1]; // d fp[0] / d p[1]
        fpd[1, 0] = d[1] * (2 * radialDistortion[0] + (4 * radialDistortion[1] + 6 * radialDistortion[2] * r2) * r2) * d[0] + 2 * tangentialDistortion[0] * d[0] + 2 * tangentialDistortion[1] * d[1]; // d fp[1] / d p[0]
        fpd[1, 1] = div + d[1] * (2 * radialDistortion[0] + (4 * radialDistortion[1] + 6 * radialDistortion[2] * r2) * r2) * d[1] + 2 * tangentialDistortion[1] * d[0] + 6 * tangentialDistortion[0] * d[1]; // d fp[0] / d p[0]

        /* Perform the Newton-Raphson step: */
        float det = fpd[0, 0] * fpd[1, 1] - fpd[0, 1] * fpd[1, 0];
        p[0] -= (fpd[1, 1] * fp[0] - fpd[0, 1] * fp[1]) / det;
        p[1] -= (fpd[0, 0] * fp[1] - fpd[1, 0] * fp[0]) / det;
    }

radialDistortion[] is the k coefficients, tangentialDistortion[] is set to 0. undistortMaxError is 1.0e-32f, undistortMaxSteps is 20.

Input:

new Vector2((mx * kinectSize.x / (meshW - 1) - depthParams.principalPoint.x) / depthParams.focalLength.x,
            (my * kinectSize.y / (meshH - 1) - depthParams.principalPoint.y) / depthParams.focalLength.y);

resulting in input points from (-0.7013, -0.5808) to (0.7013, 0.5808)

I stepped through the code at some of the bad corners, and the loop seemed to diverge instead of converge on the correct result like the other points.

Edit: turns out I made a mistake when translating the fpd section. I used radialDistortion for a few variables when I should have used tangentialDistortion. Now it works for the specific values for my Kinect, but I still get some errors in the corners with larger coefficients. I have updated the code and question to reflect this