Welcome to Fractal Forums

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
}

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);
}