// ======================================================================
// W space: Multi star - fractal sound creation tool
// (CC) 2009 Jens Koeplinger, http://www.jenskoeplinger.com/P
// License: Creative Commons Attribution "Share Alike 3.0 Unported"
// http://creativecommons.org/licenses/by-sa/3.0
// See main file "wspacemultistarsound.c" for what this means to you!
// ======================================================================
// lib-multistar-recursion.h
// ----------------------------------------------------------------------
// calculates the convergence percentage for a given point (x, y)
// in W space, with precision down to a specified threshold recursion
// depth.
//
// Example:
/*

    // calculate convergence percentage at coordinate (X, Y):

    thisPercent = squareRecursive(X, Y, 0, thresholdDepth);

      ...

    // calculate radius from (0, 0) to the outer rim of
    // thisConvPercentage at angle thisAngle, with full amplitude (1.):

    thisRadius = calculateRadius(thisConvPercentage, thisAngle, 1.);

*/
// ---------------------------------------------------------------------- 

double squareRecursive(double inReal, double inW, int inDepth, int maxDepth) {
// IN:  (inReal, inW) is the point to be squared
//      inDepth is the current recursion depth
//      maxDepth is the maximum depth
// OUT: average convergent points, [ 0. , 1. ]

    double outReal;
    double outReal2;
    double outWPlus;
    double outWMinus;
    double normPlusPlus;   // the larger of the norms on the (w) orbit result
    double normPlusMinus;  // the greater of the norms on the (w) orbit result
    double normMinusPlus;  // the greater of the norms on the (-w) orbit result
    double normMinusMinus; // the lesser of the norms on the (-w) orbit result
    double convPlus;       // 1 if (w) orbit point is convergent, 0 if divergent, 0.5 if undecided
    double convMinus;      // same for (-w) orbit point

    outReal   = inReal * inReal - inW * inW;
    outReal2  = outReal * outReal;
    outWPlus  = (2. * inReal + inW) * inW;
    outWMinus = (2. * inReal - inW) * inW;

    normPlusPlus   = sqrt( outReal2   + outReal * outWPlus    + outWPlus  * outWPlus  );
    normPlusMinus  = sqrt( outReal2   - outReal * outWPlus    + outWPlus  * outWPlus  );
    normMinusPlus  = sqrt( outReal2   + outReal * outWMinus   + outWMinus * outWMinus );
    normMinusMinus = sqrt( outReal2   - outReal * outWMinus   + outWMinus * outWMinus );

    if (inDepth == maxDepth) {
        // maximum recursion depth reached;
        // check each of the four norms and build a rough
        // estimate on behavior for infinity

        convPlus = 0.;

        if (normPlusPlus   < 1.0) convPlus = convPlus + 0.25;
        if (normPlusMinus  < 1.0) convPlus = convPlus + 0.25;
        if (normMinusPlus  < 1.0) convPlus = convPlus + 0.25;
        if (normMinusMinus < 1.0) convPlus = convPlus + 0.25;

        return convPlus;

    } else {
        // maximum recursion depth not yet reached;
        // check each point

        // 1) check (w) point

        if (normPlusMinus > 3.) {
            // must be divergent for all subsequent recursions; abort

            convPlus = 0.;

        } else if (normPlusPlus < 0.333333333333) {
            // must be convergent for all subsequent recursions; abort

            convPlus = 1.;

        } else {
            // divergence / convergence to be determined one level deeper

            convPlus = squareRecursive(outReal, outWPlus, inDepth + 1, maxDepth);

        }

        // 2) check (-w) point

        if (normMinusMinus > 3.) {
            // must be divergent for all subsequent recursions; abort

            convMinus = 0.;

        } else if (normMinusPlus < 0.333333333333) {
            // must be convergent for all subsequent recursions; abort

            convMinus = 1.;

        } else {
            // divergence / convergence to be determined one level deeper

            convMinus = squareRecursive(outReal, outWMinus, inDepth + 1, maxDepth);

        }

        // ... and return the average

        return (convPlus + convMinus) / 2.;

    }

}

// ---------------------------------------------------------------------- 

double calculateRadius(double thisConvPercentage, double thisAngle, double thisAmplitude) {
// Calculates the distance of the coorinate origin (0, 0) to the outer
// rim of the multi-star fractal's convergence percentage thisConvPercentage,
// at an angle of thisAngle (in radians).
//
// IN:  thisConvPercentage - search for outer rim at this convergence percentage
//      thisAngle - angle in radians
//      thisAmplitude - a factor of 1. or lesser, to correct for loudness
// OUT: distance (radius) from coordinate origin to outer rim of the specified
//      convergence percentage of the multi-star fractal

    double thisUnitVecX;              // real part of unit vector in direction thisAngle
    double thisUnitVecY;              // ... w part
    double prevRadiusMin;             // previous radius (minimum bound of calculation result)
    double prevRadiusMax;             // ... (maximum bound)
    double thisRadius;                // current radius to test
    double thisRadiusStep;            // absolute value of last radius correction
    double thresholdRadiusStep;       // threshold for radius corrections; once corrections in the
                                      //     radius are smaller than this threshold, abort the
                                      //     refinemenet calculation and report the result
    double thisPercentage;            // current conv percentage

    // calculation works through iterative, convergent refinement; each
    // refinement step must necessarily be smaller than the previous one (by
    // the choice of this algorithm), and always converges. Once the refinement
    // step is smaller than a certain threshold, abort the iteration and report
    // the result. The threshold is given by the resolution of the WAV file sound
    // sample (i.e. the wave table amplitude resolution, +/- 32766 integer).

    // set seed values
    thisUnitVecX   = cos(thisAngle);
    thisUnitVecY   = sin(thisAngle);
    prevRadiusMin  = 0.65;               // at the coordinate origin, we have full convergence
    prevRadiusMax  = 3.  ;               // no divergent point is further away than radius 3.
    thisRadius     = 1.  ;               // start on the unit circle
    thisRadiusStep = 1.  ;
    thresholdRadiusStep = 1. / (60000. * thisAmplitude);

    while (thisRadiusStep > thresholdRadiusStep) {
        // zoom in on the outer rim of thisConvPercentage

        thisPercentage = squareRecursive(
                             thisUnitVecX * thisRadius,
                             thisUnitVecY * thisRadius,
                             0,
                             thresholdDepth);

        if (thisConvPercentage < thisPercentage) {
            // requested convergence percentage is lower than at the current radius
            // => increase the radius (further out = lower convergence percentage)

            prevRadiusMin  = thisRadius;
            thisRadius     = (thisRadius + prevRadiusMax) / 2.;
            thisRadiusStep = thisRadius - prevRadiusMin;

        } else {
            // requested convergence is larger than convergence at the current radius
            // => decrease radius (further in = higer convergence percentage)

            prevRadiusMax  = thisRadius;
            thisRadius     = (thisRadius + prevRadiusMin) / 2.;
            thisRadiusStep = prevRadiusMax - thisRadius;

        }

    }

    return thisRadius;

}

// ----------------------------------------------------------------------