where t=0.6151732.
We used rather straightforward boundary scanning method (see e.g. D. Saupe, Efficient computation of Julia sets and their fractal dimension, Physica D, v.28 n.3, p.358-370, Oct. 1987 [1]). C++ source code follows. It's available under CC BY-SA, GFDL and GPLv2+.
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <iostream>
#include <complex>
inline double max(double a,double b)
{
return a>b?a:b;
}
using namespace std;
const double t=0.6151732;
const complex<double> rot=exp(2*M_PI*complex<double>(0,1)*t);//e^{2 \pi i t}
complex<double> f(complex<double> z)
{
// f(z) = e^{2 \pi i t}\cdot \frac{z^2(z - 4)}{1 - 4z},
return rot*(z*z*(z-4.))/(complex<double>(1,0)-4.*z);
}
int main()
{
//z=x+iy;
double xmin=-4.0,xmax=8.0,deltax,ymin=-3,ymax=3,deltay;
int precx=400, precy=200;
deltax=(xmax-xmin)/precx;
deltay=(ymax-ymin)/precy;
int iterations=4000;
double threshold=iterations/30*sqrt(deltax*deltax+deltay*deltay);
double x,y;
complex<double> zlt,zlb,zrt,zrb;
int i,k,l;
double maxdistance;
int color, MaxColor=15;
double infty=1.E4;
bool inf_zlt,inf_zlb,inf_zrt,inf_zrb;
cout << "P2" << endl << precx << " " << precy << endl << MaxColor << endl;
for(k=0;k<precy;k++)
{
y=ymax-k*deltay;
cerr << y << endl;
for(l=0;l<precx;l++)
{
x=xmin+l*deltax;
zlt=complex<double>(x-deltax/2,y+deltay/2);
zlb=complex<double>(x-deltax/2,y-deltay/2);
zrt=complex<double>(x+deltax/2,y+deltay/2);
zrb=complex<double>(x+deltax/2,y-deltay/2);
inf_zlt=false;
inf_zlb=false;
inf_zrt=false;
inf_zrb=false;
for(i=0;i<iterations;i++)
{
if(!inf_zlt){zlt=f(zlt);}
if(!inf_zlb){zlb=f(zlb);}
if(!inf_zrt){zrt=f(zrt);}
if(!inf_zrb){zrb=f(zrb);}
if(abs(zlt)>infty)
{
zlt=infty;
inf_zlt=true;
}
if(abs(zrt)>infty)
{
zrt=infty;
inf_zrt=true;
}
if(abs(zlb)>infty)
{
zlb=infty;
inf_zlb=true;
}
if(abs(zrb)>infty)
{
zrb=infty;
inf_zrb=true;
}
}
maxdistance=max(abs(zlt-zlb),max(abs(zlt-zrt),max(abs(zlt-zrb),max(abs(zlb-zrt),max(abs(zlb-zrb),abs(zrt-zrb))))));
if(maxdistance>threshold*8)
{
color=0;
}
else
{
if(maxdistance<threshold/8.)
{
color=MaxColor;
}
else
{
color=MaxColor-(int(log(maxdistance/threshold)/log(8)*7.)+7);
}
}
//if (x>0 && y>0) color=MaxColor-color; /* check the orientation of Z-plane by marking first quadrant */
cout << color << " ";
}
cout << endl;
}
return 0;
}