#include <stdio.h>
#include <stdlib.h>
#include <math.h>

static double fact( int n ) {  // The factorial function!
	double f = 1;
	while(n>1)
		f *= n--;
	return f;
}

static inline unsigned char f2b( double f ) {  // float to byte
	if(f<0) f = 0; if(f>1) f = 1;
	int i = f*256; if(i>255) i = 255;
	return i;
}

void HSV2RGB(double h, double s, double v, unsigned char rgb[3] ) {  // the classic hue scale
	if (s == 0) {
		rgb[0] = rgb[1] = rgb[2] = f2b(v);
	} else {
		double v_h = h * 6;
		double v_i = floor(v_h);
		double v_1 = v * (1 - s);
		double v_2 = v * (1 - s * (v_h - v_i));
		double v_3 = v * (1 - s * (1 - (v_h - v_i)));
		double v_r,v_g,v_b;
		     if (v_i == 0) {v_r = v;   v_g = v_3; v_b = v_1;}
		else if (v_i == 1) {v_r = v_2; v_g = v;   v_b = v_1;}
		else if (v_i == 2) {v_r = v_1; v_g = v;   v_b = v_3;}
		else if (v_i == 3) {v_r = v_1; v_g = v_2; v_b = v  ;}
		else if (v_i == 4) {v_r = v_3; v_g = v_1; v_b = v  ;}
		else               {v_r = v;   v_g = v_1; v_b = v_2;};
		
		rgb[0] = f2b(v_r);
		rgb[1] = f2b(v_g);
		rgb[2] = f2b(v_b);
	}
}

double zernike_polynomials( int m, int n, double ro, double th ) {  // the polynomial hitself
	double Rmnro = 0;
	bool even = m>=0;
	if(m<0) m = -m;
	if( (n-m)%2 ) return 0;

	for(int k=0;k<=(n-m)/2;++k) {
		Rmnro += pow(ro,n-2*k)*
			     ( pow(-1.0,k) * fact(n-k) ) /
			     ( fact(        k) *
			       fact((n+m)/2-k) * 
			       fact((n-m)/2-k)
			     );
	}

	if(even) return Rmnro * cos(m*th);
	else     return Rmnro * sin(m*th);
}

void main() {
	const double NO_VALUE = 42; const int B  = 64;
	const int SX = 1024;
	const int SY = 1024;
	
	const int LX = 3; const int LY = 7;
	int pairs[21][2] =
	{
		{ 0,0},
		{-1,1},{ 1,1},
		{-2,2},{ 0,2},{ 2,2},
		{-3,3},{-1,3},{ 1,3},{ 3,3},
		{-4,4},{-2,4},{ 0,4},{ 2,4},{ 4,4},
		{-5,5},{-3,5},{-1,5},{ 1,5},{ 3,5},{ 5,5},
	};

	unsigned char * img = new unsigned char[LX*SX*LY*SY*3];

	int ix,iy,q = 0;
	for(iy=0;iy<LY;++iy) for(ix=0;ix<LX;++ix) {
		int m = pairs[q][0];
		int n = pairs[q][1];
		for(int j=0;j<SY;++j) {
			double y = double(SY/2-j)/(SY/2-B);
			for(int i=0;i<SX;++i) {
			   double x = double(SX/2-i)/(SX/2-B);
			   double ro = sqrt(x*x+y*y);
			   double th = atan2(y,x);	
			   double z = NO_VALUE;
			   if(ro<=1)
					z = zernike_polynomials(m,n,ro,th);
			   unsigned char * rgb = img+3*( i+SX*(ix+LX*(j+SY*iy)) );
				if(z == NO_VALUE )
					rgb[0] = rgb[1] = rgb[2] = 255;
				else {
					if(z>=0) z =  pow( z,0.8);	// little enance
					else     z = -pow(-z,0.8);
					z = (z+1)/2;			// normalize
					if(z<0) z = 0;
					if(z>1) z = 1;
					
					HSV2RGB( 2.0*z/3.0, 0.8, 0.9, rgb );
				}
			}
		}
		printf("%03d%%\n",q*100/(LX*LY)); ++q;
	}

	FILE * fp = fopen("c:\\temp\\zernike.ppm","wb");
	fprintf(fp,"P6\n%d %d\n255\n",SX*LX,SY*LY);
	fwrite(img,1,LX*SX*LY*SY*3,fp);
	fclose(fp);
	delete[] img;
}