END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

// Mandelbrot iteration of a single point
// tries to compute angle of external ray passing through that point (buggy!)
unsigned iterate(const double cx, const double cy, const long maxiter) {
    unsigned i = 0;                 // iteration count
    double zx = 0.0, zy = 0.0;      // complex z
    double z2x = 0.0, z2y = 0.0;    // z real and imaginary parts squared
    
    unsigned angle = 0;             // accumulates fractional angle
    unsigned curBit = 0x80000000UL; // current bit of binary fraction
    
    do {
        i++;
        if (i > maxiter) {
            return -1;      // alternative loop end
        }
        
        z2x = zx*zx;        // squared components
        z2y = zy*zy;

        zy = 2.0*zx*zy;     // complex square, in place
        zx = z2x - z2y;
        
        zx = zx + cx;       // plus constant
        zy = zy + cy;
        
        // accumulate another bit of fractional angle
        if (zy < 0.0) {
            angle |= curBit;
        }
        curBit = curBit >> 1;   // next bit
        
    } while ((zx*zx + zy*zy) < 10000.0);
    
    // keep refining the angle for escaped points
    while (curBit > 0) {
        z2x = zx*zx;        // squared components
        z2y = zy*zy;
        
        zy = 2.0*zy*zx;     // complex square
        zx = z2x - z2y;
        
        zx += cx;           // plus constant
        zy += cy;
        
        // accumulate another bit of fractional angle
        if (zy < 0.0) {
            angle |= curBit;
        }
        curBit = curBit >> 1;   // next bit
    }
    
    return angle;
}

MMF3-FieldLines {
;
; This version of calculating field lines will work reasonably well with both
; Mandelbrot and Julia Sets for divergent fractals with a divergence that is
; close to being a positive integer >=2 but can also produce quite interesting
; results with fractals that do not fit that criteria - just not rendering the
; field lines correctly in other cases :-)
;
; Thanks to Chris Hayton for showing me how to fix the "off-by-one" problems.
;
global:
  float twopi = 2.0*#pi
  float dtwopi = 0.5/#pi
  float dp = 1.0/@power
  float m = @power*dtwopi
  float dp2 = 1.0/(2.0+@power) ; Changed from originally trying 1/power
init:
  complex zv[#maxiter]
  int i = 0
  float a = 0.0
  float b = 0.0
  float c = 0.0
  float s = 0.0
  if @julia
    zv[i] = #z
    i = i + 1
  endif
loop:
  zv[i] = #z
  i = i + 1
final:
  if (s = atan2(#z))<0
    s = s + twopi
  endif
  while i > @skipiter
    i = i - 1
    if @version==0
      b = (@power*atan2(zv[i]))%twopi
    else
      b = atan2(zv[i]^@power)
    endif
    if b < 0.0
      b = b + twopi
    endif
    c = c + b - s
    if c >= #pi
      s = s + twopi
      if @fixit || @accident
        c = c - twopi
      endif
    elseif c < -#pi
      s = s - twopi
      if @fixit || @accident
        c = c + twopi
      endif
    endif
    if @fixit
      c = dp2*c
    elseif !@accident
      c = 0.0
    endif
    if (a = atan2(zv[i])) < 0.0
      a = a + twopi
    endif
    s = dp*(s + twopi*floor(a*m))
  endwhile
  s = s*dtwopi
  if @fixskip
    while s<0.0
      s = s + 1.0
    endwhile
    while s>=1.0
      s = s - 1.0
    endwhile
  endif
  #index = s
default:
  title = "MMF3 Field Lines (use high bailouts)"
  heading
    caption = "Information"
;    text = "This colouring is a little touchy around certain values, if you \
;            get obvious lines or dots that shouldn't be there but the 'Degree \
;            of Divergence' is set properly for the main formula and you are \
;            using a high bailout of say 1e20 or higher then try \
;            adjusting the location very slightly. For example if vertical \
;            lines appear on an unrotated fractal then try changing the x \
;            (real) location very slightly, if horizontal lines appear then try \
;            changing the y (imag) location very slightly."
  endheading
  heading
  endheading
  param version
    caption = "Version"
    enum = "Original" "Alternative"
    default = 1
    hint = "Allows you to modify one of the calculations in the formula that \
            is particularly sensitive to small errors. Ideally both methods \
            should produce identical results but they don't, so if you get \
            visible errors using one method then try the other (assuming you \
            have the divergence set correctly and are using a high bailout)."
  endparam
  param power
    caption = "Degree of divergence"
    default = 2.0
    hint = "This is an estimate of the divergence of the main formula, e.g. 2 \
            for z^2+c, 3 for z^3+c etc."
  endparam
  param julia
    caption = "Julia ?"
    default = false
    hint = "Enable for correct rendering of Julia Sets."
  endparam
  param skipiter
    caption = "'Start' iteration"
    default = 0
    min = 0
    hint = "You can use this to attempt to get better detail at higher \
            iteration depths, note that it only works as intended with 'Color \
            Density' as a whole integer >=1."
  endparam
  param fixskip
    caption = "Fix 'Start' iteration"
    default = false
    hint = "Fixes a problem when the start iteration is non-zero and maybe in \
            other cases - specifically if you get solid areas of palette \
            colour zero when the Transfer Function is Linear this may fix the \
            problem."
  endparam
  param fixit
    caption = "Fix Field Lines"
    default = false
    hint = "Enable this to make the formula much more accurate at rendering \
            the field lines correctly. It's a fudge produced after trial and \
            error but is quite effective at correcting some areas."
  endparam
  param accident
    caption = "Happy Accident"
    default = false
    hint = "You may find the colouring useful when this is enabled. This \
            is available as I found the results when initially adding the \
            'Fix it' option quite interesting."
    visible = !@fixit
  endparam
}

void getGradient(long double cx, long double cy, long double &xDer,
                 long double &yDer, long double &gzx, long double &gzy,
                 long double &length, unsigned int maxIter) {
        // based on gradient equation found at:
        // http://linas.org/art-gallery/escape/ray.html
    long double zx = 0.0L;
    long double zy = 0.0L;
    xDer = 0.0L; yDer = 0.0L; // initialize values
    long double xTemp;
    unsigned int n = 0;
    long double huge = std::numeric_limits<float>::max();
    while ((zx*zx + zy*zy) < huge && n < maxIter) {
            // compute new derivative:
        xTemp = 2.0L*(zx*xDer - zy*yDer) + 1.0L;
        yDer = 2.0L*(zx*yDer + zy*xDer);
        xDer = xTemp;
            // compute new z value:
        xTemp = zx*zx - zy*zy + cx;
        zy = 2.0L*zx*zy + cy;
        zx = xTemp;
        n++;
    }
    length = std::sqrt(zx*zx + zy*zy);
    long double m = n + 1 - (std::log(std::log(length)) / std::log(2.0L));
    long double f = pow(2.0L,-m); // Douady-Hubbard potential
    long double Dx = zx*(xDer/2.0) - zy*(-yDer/2.0); // real part of zn * Dn
    long double Dy = zx*(-yDer/2.0) + zy*(xDer/2.0); // imag part of zn * Dn
    gzx = f*Dx / (length*length*std::log(length)); // x-gradient for cx,cy
    gzy = f*Dy / (length*length*std::log(length)); // y-gradient for cx,cy
}

double externalAngle(long double cx, long double cy, std::size_t maxIter) {
        // uses 4th order Runge-Kutta integration to find external angle
    long double huge = 640000.0; // bigger is better, but also slower
    long double ckx = cx;
    long double cky = cy;
    long double ckxTemp, ckyTemp;
    long double xDer,yDer,length,derLength,dist,gzx,gzy,gzxA,gzyA,gzxB,gzyB,gzxC,gzyC,gzxD,gzyD,gzlength;
    unsigned int k = 0;
    unsigned int maxStep = maxIter*2; // make sure we try to walk far away
    while((ckx*ckx + cky*cky) < huge && k++ < maxStep) {
        // integrate with Runge-Kutta outwards along gradient:
        getGradient(ckx,cky,xDer,yDer,gzxA,gzyA,length,maxIter);
        derLength = std::sqrt(xDer*xDer + yDer*yDer);
        dist = 0.5L * std::log(length) * length / derLength;
        dist = std::min(16.0L,0.5*dist); // reduce step size
        gzlength = std::sqrt(gzxA*gzxA + gzyA*gzyA);
        ckxTemp = ckx + (0.5 * dist * gzxA/gzlength); // walk to midpoint
        ckyTemp = cky + (0.5 * dist * gzyA/gzlength); // using A
        getGradient(ckxTemp,ckyTemp,xDer,yDer,gzxB,gzyB,length,maxIter);
        gzlength = std::sqrt(gzxB*gzxB + gzyB*gzyB);
        ckxTemp = ckx + (0.5 * dist * gzxB/gzlength); // walk to midpoint
        ckyTemp = cky + (0.5 * dist * gzyB/gzlength); // using B
        getGradient(ckxTemp,ckyTemp,xDer,yDer,gzxC,gzyC,length,maxIter);
        gzlength = std::sqrt(gzxC*gzxC + gzyC*gzyC);
        ckxTemp = ckx + (dist * gzxC/gzlength); // walk full step
        ckyTemp = cky + (dist * gzyC/gzlength); // using C
        getGradient(ckxTemp,ckyTemp,xDer,yDer,gzxD,gzyD,length,maxIter);
        gzx = (1.0L/6.0L) * (gzxA + 2.0L*(gzxB + gzxC) + gzxD); // average the values together
        gzy = (1.0L/6.0L) * (gzyA + 2.0L*(gzyB + gzyC) + gzyD);
        gzlength = std::sqrt(gzx*gzx + gzy*gzy); // get length to normalize
        ckx += dist * gzx/gzlength; // normalize gradient and scale by dist
        cky += dist * gzy/gzlength;
    }
    return (std::atan2(cky,ckx));
}

dist = std::min(16.0L,0.5*dist);

dist = std::min(1.0L,0.05*dist);