Need help with distance estimation formula for the Phoenix fractal

[top-quark]

I'm writing a fractal library in JavaScript for an online fractal exploration application (not yet live - I'm writing the documentation).

One of the formulae I'm offering is the Phoenix fractal. It's not looking too bad:


phoenx+0.5667-0.5P by c.c.williams, on Flickr

However, if I turn off exterior colouring and enable boundary tracing, things don't look so good:


phoenx-distance-estimate by c.c.williams, on Flickr

It looks wrong (doesn't exactly scream "nowhere dense") but not completely wrong (the shape is correct).

Formula: z_k+1 = z_k * z_k + c + P * z_k-1

Here's my algorithm:

Parameters: z, c, P

z_ (im(z), re(z)),    // side-by-side view
zlast = 0,
iter = 0,
eV = 256,
dz = 1,               // Cumulative differential of f(z_)
dzp = 0,              // Cumulative differential of P * f(zlast)
rsquared = re(z_) * re(z_) + im(z_) * im(z_);

while (iter < MAX_ITER && rsquared < eV) {
    _z = z_ * z_ + c + P * zlast;
    dz = 2 * z_ * dz + dzp;
    dzp = 2 * P * zlast * dzp;
    zlast = z_;
    z_ = _z;
    rsquared = re(z_) * re(z_) + im(z_) * im(z_);
}

if (iter < MAX_ITER) {
    r = |z_|;
    d = |dz];
    distance = log(r * r) * r / d;
    if (TRACE_BOUNDARY AND distance < PIXEL_WIDTH / 20) {
        return BLACK;
    }
    else if (EXTERIOR_COLOURING) {
         return COLOUR;
    }
}
return WHITE;



If any knowledgeable person can point me in the direction of where I'm going wrong in incorporating P * zlast into the cumulative differential of ∂f(z) / dz, I'd be very grateful.

(P.S. Distance estimate works fine for Mandelbrot / Multibrot / Burning Ship and corresponding Julia sets).

[knighty]

Formula: z_k+1 = z_k * z_k + c + P * z_k-1

Are you sure this is the formula for the phoenix fractal? isn't it:
   z_k+1 = z_k * z_k + Re(c) + Im(c) * z_k-1 ?
(Edit: unless c=Re(c) and P=Im(c) of course  )
This transformation being not conformal here is a routine that uses the jacobian matrix instead to compute DE for julia sets:

Code:

phoenix(x,y,cx,cy,mi){
   //the computation of the DE is for julia sets here-> cx and cy are constant
   lx=0; ly=0;//z[n-1]
   dlxx=0; dlxy=0; dlyx=0; dlyy=0;//Dz[n-1] (jacobian matrix)
   dxx=1; dxy=0; dyx=0; dyy=1;//Dz[n]
   r2=x*x+y*y;
   for(i=0; i<mi && r2<1024; i++){
      //Compute Dz[n+1]
	  ndxx= 2*(x*dxx-y*dyx)+cy*dlxx;
      ndxy= 2*(x*dxy-y*dyy)+cy*dlxy;
      ndyx= 2*(x*dyx+y*dxx)+cy*dlyx;
      ndyy= 2*(x*dyy+y*dxy)+cy*dlyy;
      //set Dz[n-1]:=Dz[n]
	  dlxx=dxx; dlxy=dxy; dlyx=dyx; dlyy=dyy;
	  //set Dz[n]:=Dz[n+1]
      dxx=ndxx; dxy=ndxy; dyx=ndyx; dyy=ndyy;
      //compute z[n+1]
	  nx= x*x-y*y+cx+cy*lx;
      ny= 2*x*y+cy*ly;
      //set z[n-1]:=z[n]
	  lx=x;ly=y;
	  //set z[n]:=z[n+1]
      x=nx; y=ny;
	  //Modulus squared of z[n]
      r2=x*x+y*y;
   }
   r=sqrt(r2);//modulus of z[n]
   dr=max(sqrt(dxx*dxx+dxy*dxy),sqrt(dyx*dyx+dyy*dyy));//a norm of Dz[n]
   return max(0,0.5*r*log(r)/dr);//Distance estimate
}

[knighty]

Oops! just realized that the phoenix julia set transformation is conformal just because c Is constant.

Code:

phoenix(x,y,cx,cy,mi){
   //x=cx; y=cy;
   lx=0; ly=0;
   dlx=0; dly=0;
   dx=1; dy=0;
   r2=x*x+y*y;
   for(i=0; i<mi && r2<10240; i++){
      ndx= 2*(x*dx-y*dy)+cy*dlx;
      ndy= 2*(x*dy+y*dx)+cy*dly;
      dlx=dx; dly=dy;
      dx=ndx; dy=ndy;
      nx= x*x-y*y+cx+cy*lx;
      ny= 2*x*y+cy*ly;
      lx=x;ly=y;
      x=nx; y=ny;
      r2=x*x+y*y;
   }
   r=sqrt(r2);
   
   dr=sqrt(dx^2+dy^2);
   max(0,0.5*r*log(r)/dr);
}

[top-quark]

Thanks for the response. The formula I'm using uses two complex constants (Ushiki) rather than than real constants (Stevens). When P is 0, you have the quadratic Julia formula z -> z² + c. However, since an imaginary component in either c or P results in crap as far as I can tell, I may well change it to use real constants. I will
certainly try the dz calculation that you were kind enough to provide.

The examples, as you probably guessed, are using c=0.5667 and P=-0.5.

[top-quark]

Works absolutely perfectly. Muchas gracias!

[cKleinhuis]

so, what about a 3d version then ?

[top-quark]

I'll have to teach myself WebGL first, so maybe a project for the future. A decent JavaScript engine is impressively fast at floating point calculations and web workers suppress slow script warnings so it should be doable.

Right now, I'm trying to get boundary tracing for Newtonian fractals working.