File:Herman-ring-1.png

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