Welcome to Fractal Forums

vector<Element>* Newton::calculateNoOrbit(vector<Point>* points)
{
    vector<Point>::const_iterator i = points->begin();

    int iterations;
    double zx;
    double zy;
    double z2x;
    double z2y;
    double z3x;
    double z3y;
    double zx_old = DBL_MAX;
    double zy_old = DBL_MAX;
    double zx_older = DBL_MAX;
    double zy_older = DBL_MAX;
    double zx_delta;
    double zy_delta;
    double denominator;

    double radius;
    double modulus;
    double phi;

    const double exponent2 = exponent - 1.0;
    const double temp1 = 0.5 * exponent;
    const double temp2 = 0.5 * exponent2;

    getCoord(i->x, i->y, zx, zy);

    for(iterations = 0; iterations < maxIterations; ++iterations)
    {
        // Calculate |z - z_old| for bailout check
        zx_delta = zx - zx_old;
        zy_delta = zy - zy_old;

        zx_delta *= zx_delta;
        zy_delta *= zy_delta;

        // Check |z - z_old| against bailout
        if(zx_delta + zy_delta < bailout)
        { break; }

        // Update z_older and z_old
        zx_older = zx_old;
        zy_older = zy_old;
        zx_old = zx;
        zy_old = zy;

        // Calculate z^exponent - 1
        radius = zx * zx + zy * zy;
        modulus = exp(log(radius) * temp1);
        phi = exponent * atan2(zy, zx);

        z3x = modulus * cos(phi) - 1.0;
        z3y = modulus * sin(phi);

        // Calculate exponent * z^(exponent - 1)
        modulus = exp(log(radius) * temp2);
        phi = exponent2 * atan2(zy_old, zx_old);

        z2x = modulus * cos(phi) * exponent;
        z2y = modulus * sin(phi) * exponent;

        // Calculate the division and subtract from z
        denominator = z2x * z2x + z2y * z2y;

        if(denominator == 0.0)
        { break; }

        zx -= (z3x * z2x + z3y * z2y) / denominator;
        zy -= (z3y * z2x - z3x * z2y) / denominator;
    }

    results.clear();

    // Calculate smooth color
    zx_delta = zx - zx_old;
    zy_delta = zy - zy_old;
    double m1 = log(sqrt(zx_delta * zx_delta + zy_delta * zy_delta));

    zx_delta = zx_old - zx_older;
    zy_delta = zy_old - zy_older;
    double m2 = log(sqrt(zx_delta * zx_delta + zy_delta * zy_delta));

    double fractional = iterations + (log(-0.5 * log(bailout)) - log(-0.5 * m1)) / log( m1 / m2 );
    writeResult(i->x, i->y, iterations, fractional);

    return &results;
}