I also remember a thread here in fractalforums about that fractal. Unfortunately, I can't find it by now.
See also:
LINK describing an escape time definition.
I've also written a DE for it in the hope that it's 3D generalisation would give good looking fractal... That was not the case.
//Use with evaldraw (http://advsys.net/ken/download.htm)
(x,y){
static binit=1;
static frm=-1,mx,my,mr2;
if(frm<numframes){
frm=numframes;
mx=(mousx-0.5*xres)*0.005;
my=(0.5*yres-mousy)*0.005;
mr2=deifs(mx,my);mr2*=mr2;
}
if((x-mx)^2+(y-my)^2<mr2) return 0;
deifs(x,y)^0.25
}
static maxiter=100,maxdepth=30;
static tcx=0,tcy=1,tcr=sqrt(1);//tcr=sqrt(tcx^2+tcy^2)//is this thing invariant, the same for any tcx and tcy?
enum {Npq=256};
deifs(x,y){
r2=x*x+y*y;inverted=0;
if(r2<1){r2=1/r2;x*=r2;y*=r2;inverted=1;}
if(r2<1.75) return 0.04;
trx=x;
try=y;
s=sqrt(r2);r2=s;
bvr=tcr/(s-1);//bvr=min(bvr,3);
tis=1/s;
//priority queue... the most stupid possible :o))
auto pqx[Npq];//x
auto pqy[Npq];//y
auto pqs[Npq];//1/accumulated_scale
auto pqd[Npq];//distance
auto pqh[Npq];//depth
stptr=0;//pointer in the priority queue
bvrm=100;//controls the smoothness of the DE field (lower values faster, higher values smoother)
{//first instances
x0=tcx;y0=tcy;//initial value, faster than beginning at (0,0)
pqx[stptr]=x0; pqy[stptr]=y0;
pqs[stptr]=tis;//instead of 1 because of the initial value of (1,0)
pqd[stptr]=tcr-bvr;
pqh[stptr]=0;
stptr++;
}
bNoExit=1;i=0;ii=0;
while(bNoExit && i<maxiter){
//extract nearest element
{
depth=pqh[stptr-1];
if(depth>maxdepth) break;
stptr--;
x=pqx[stptr]; y=pqy[stptr];s=pqs[stptr];ad=pqd[stptr];
}
//generate and insert its children
nx=x*trx-y*try;
ny=x*try+y*trx;
nc=-1;
for(j=0;j<2;j++){
//generate child j
pqx[stptr]=nx+nc*tcx;
pqy[stptr]=ny+nc*tcy;
pqs[stptr]=s*tis;
pqd[stptr]=(sqrt((pqx[stptr])^2+(pqy[stptr])^2)-bvr)*pqs[stptr];
pqh[stptr]=depth+1;
nc=1;
//insert it at its place
ptr=stptr;
imm=0;
while(ptr>0){
if(pqd[ptr]>pqd[ptr-1]){//swap
k=pqx[ptr]; pqx[ptr]=pqx[ptr-1]; pqx[ptr-1]=k;
k=pqy[ptr]; pqy[ptr]=pqy[ptr-1]; pqy[ptr-1]=k;
k=pqs[ptr]; pqs[ptr]=pqs[ptr-1]; pqs[ptr-1]=k;
k=pqd[ptr]; pqd[ptr]=pqd[ptr-1]; pqd[ptr-1]=k;
k=pqh[ptr]; pqh[ptr]=pqh[ptr-1]; pqh[ptr-1]=k;
ptr--;imm++;
} else break;
}
stptr++;
}
d=pqd[stptr-1];x=pqx[stptr-1];y=pqy[stptr-1];s=pqs[stptr-1];
bNoExit=(x*x+y*y<bvrm*bvr*bvr);
i++;
}
DE=0.5*pqd[stptr-1]/bvr;
if(inverted)
return DE;//good for the interior!
return r2^2/(1/DE+r2);//deduced from the interior distance
}