Function complex_iterative_inverse_impl

src/transformations/horner.cpp:440–490

}

static PJ_UV complex_iterative_inverse_impl(PJ *P, const HORNER *transformation,
                                            PJ_UV position) {

    double n, e;
    // in this case fwd_origin and any existing flipping needs to be added in
    // the end
    e = position.u;
    n = position.v;

    if (coords_out_of_range(P, transformation, n, e)) {
        return generate_error_coords();
    }

    {
        // complex real part corresponds to Northing, imag part to Easting
        const double tol = transformation->inverse_tolerance;
        const std::complex<double> dZ(n - transformation->fwd_c[0],
                                      e - transformation->fwd_c[1]);
        std::complex<double> w0(0.0, 0.0);
        int loops = 32; // usually converges really fast (1-2 loops)
        bool converged = false;
        while (loops-- > 0 && !converged) {
            // sum coefficient pointers from back to front until the first
            // complex pair (fwd_c0+i*fwd_c1)
            const double *c = transformation->fwd_c;
            PJ_UV en = {w0.imag(), w0.real()};
            en = complex_horner_eval(transformation->order, c, en, 1);
            std::complex<double> det(en.v, en.u);
            std::complex<double> w1 = dZ / det;
            converged = (fabs(w1.real() - w0.real()) < tol) &&
                        (fabs(w1.imag() - w0.imag()) < tol);
            w0 = w1;
        }
        // if loops have been exhausted and we have not converged yet,
        // we are never going to converge
        if (!converged) {
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM);
            position = generate_error_coords();
        } else {
            double E = w0.imag();
            double N = w0.real();
            if (transformation->uneg)
                E = -E;
            if (transformation->vneg)
                N = -N;
            position.u = E + transformation->fwd_origin->u;
            position.v = N + transformation->fwd_origin->v;
        }
        return position;
    }
}

static void complex_horner_forward_4d(PJ_COORD &point, PJ *P) {
    const HORNER *transformation = reinterpret_cast<const HORNER *>(P->opaque);