#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "bitmaplib.h"

/*
	Henon phase attractor, simple example
*/

#define NX 2000
#define NY 2000
long N = 1000000L;
#define SCALE (NX / 3)

typedef struct {
	double x,y;
} XY;

int main(int argc,char **argv)
{
	long n,nnew=0;
	int ix,iy,counter=0;
	XY p0={0,0},p1;
	double angle,ca,sa,r=0;
	BITMAP4 *image,white={255,255,255,0},black={0,0,0,0};
	char fname[64];
	FILE *fptr;

	// Provide the parameter a, in radians
	if (argc < 2) {
		fprintf(stderr,"Usage: %s angle\n",argv[0]);
		exit(-1);
	}
	angle = atof(argv[1]);

	// Create the image
	image = Create_Bitmap(NX,NY);
	Erase_Bitmap(image,NX,NY,white);

	// Sin and cosine terms
	ca = cos(angle);
	sa = sin(angle);

	for (n=0;n<=N;n++) {
		// Reset the seed every now and then, range +-1 on each axis
		// Or if the point has gone to infinity
		r = p0.x*p0.x + p0.y*p0.y;
		if (r > 1e32 || n % (N/1000) == 0) {
         p0.x = 2*drand48() - 1;
         p0.y = 2*drand48() - 1;
			nnew = 0;
		}

		// Iterate
		p1.x = p0.x * ca - (p0.y - p0.x*p0.x) * sa;
		p1.y = p0.x * sa + (p0.y - p0.x*p0.x) * ca;
		if (nnew > 10) { // Found the attractor
			ix = p1.x * SCALE + NX/2; // Mapping to pixel coordinates
			iy = p1.y * SCALE + NY/2;
      	Draw_Pixel(image,NX,NY,ix,iy,black);
		}
		p0 = p1;

		// Write every 10th image to a file
		if (n > 0 && n % (N/10) == 0) {
			sprintf(fname,"henon_%02d.tga",counter);
			if ((fptr = fopen(fname,"wb")) == NULL) {
				fprintf(stderr,"Unable to open image file\n");
				exit(0);
			}
			Write_Bitmap(fptr,image,NX,NY,12);
			fclose(fptr);
			counter++;
		}

		nnew++;
	}

	exit(0);
}